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( ^ When a film of a liquid suspension of nanoparticles or a polymer solution is deposited on a surface, it may dewet from the 
surface and as the solvent evaporates the solute particles/polymer can be deposited on the surface in regular line patterns. In this 
tkQ paper we explore a hydrodynamic model for the process that is based on a long-wave approximation that predicts the deposition 
of irregular and regular line patterns. This is due to a self-organised pinning-depinning cycle that resembles a stick-slip motion 

"^Lh of the contact line. We present a detailed analysis of how the line pattern properties depend on quantities such as the evaporation 
rate, the solute concentration, the Peclet number, the chemical potential of the ambient vapour, the disjoining pressure, and the 

CO intrinsic viscosity. The results are related to several experiments and to depinning transitions in other soft matter systems. 
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1 Introduction 

Many production processes employed in the chemical, phar- 
maceutical and other industries, involve a wet phase where 
films or drops of a solution or suspension are applied to a solid 
or liquid surface with the aim of producing a homogeneous or 
structured layer of the solute on the surface. Prominent exam- 
ples are printing, painting and coating processes where a wet 
ink or paint is used. The carrier fluid (solvent) then evaporates 
and leaves all the originally dissolved non-volatile material be- 
hind in various deposition patterns on the substrate. 

If the aim of these processes is to produce a large-scale 
homogeneous layer, ideally one would instantaneously pro- 
duce an extended stable homogenous film of solution that then 
evaporates slowly and homogeneously. Although, this may be 
achieved to some extent, e.g. by spin coating, this situation is 
rather exceptional. In practice, what generally happens is that 
even while the wet film is still being deposited onto one part 
of the substrate, it is already dry on other parts. Many tech- 
niques exist for applying such films, that are relevant for dif- 
ferent surface types and surface geometries, such as painting, 
blade coating, deposition of individual drops, spray coating, 
and processes involving liquid menisci. 

Historically, the main aim of these processes was to produce 
homogeneous coatings, so defects like holes or patterned areas 
would not have been desirable. However, this has changed as 
nowadays wet deposition may also be used to produce pat- 
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terned functional layers via the self-organisation of the so- 
lute. Although a large variety of different apparatus geome- 
tries exists, the basic physical question always is What types 
of deposition patterns exist and how can the self-organisation 
processes involved informing the patterns be controlled? A 
close inspection of the individual systems shows that in most 
of them, the formation of deposition patterns occurs in a lo- 
calised zone at receding contact lines that move under the cou- 
pled influence of solvent evaporation, convective motion of the 
solution, capillarity and wettability. 

In this paper we present a long-wave hydrodynamic thin 
film model that focuses on such a moving contact line. Some 
of the main features of our model were briefly presented by 
Frastia et alM*. However, before we discuss the model and the 
wide range of detailed results that we have obtained, we first 
briefly review some of the relevant experimental systems and 
theoretical approaches that exist in the literature. 

Interest in deposition patterns has markedly increased over 
the last decade, since Deegan and co-workers' detailed inves- 
tigations of the "coffee-stain effect", i.e., the solute deposition 
patterns that are left behind by a receding three-phase contact 
line of an evaporating drop of a suspension upon a smooth 
solid substrate.^™ In particular Deegan 3 describes and analy- 
ses a wide range of patterns: cellular and lamellar structures, 
single and multiple rings, and Sierpinski gaskets. The creation 
of multiple rings through a stick-slip front motion of colloidal 
liquids was also observed. 5 6 Since then, interest widened and 
now encompasses all phenomena that accompany evapora- 
tive and convective dewetting of colloidal suspensions and 
polymer or macromolecular solutions in a number of differ- 
ent geometries. Other early examples are the investigation 
by Parisse and Allain of the shape changes that drops of col- 
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loidal suspension undergo when they dry and the creation 
of semiconductor nanoparticle rings through a similar deposi- 
tion process.! 9 -! Other observed drying structures includ e cra ck 
patterns, 10 chevron patterns,^ and branched patterns JEHU i n 
fact, crack patterns in sol-gel processes had already been stud- 
ied somewhat earlierJSEEl however, we do not consider them 
here. The related concept of using evaporation at contact lines 
to assemble colloidal particles or proteins into crystals actu- 
ally has a longer history - see for example the discussions and 
reviews by Davies et alJ^, Denkov et al. 17 , Adachi and Na- 
gayama^H Maenosono et alJ^, Kinge et al.&H. 

Generally, the evaporation of a macroscopic drop does not 
usually result in the deposition of regular concentric ring pat- 
terns that could potentially be employed to fabricate devices, 
but rather results in irregular patterns of rugged rings and 
lines (see, e.g., Deegan^, Adachi et al. 5 ). This changes 
when the experiments are performed in a controlled way on 
smaller scales: Recently, both polymer solutions 21 " 23 ! and 
(nano)particle suspensionsl^SEH have been employed in var- 
ious small-scale geometries where one is able to exercise 
greater control over the contact line as it recedes due to evapo- 
ration. As a result, strikingly regular line patterns are created, 
where the deposited structures show typical distances rang- 
ing from 10-100/im. Line patterns can be parallel or per- 
pendicular to the receding contact line and are produced in 
a robust repeatable manner in extended regions of parameter 
space. Besides the lines, a variety of other patterns may also be 
found, including undulated stripes, interconnected stripes;^ 
ladder structures, i.e. superpositions of perpendicular and par- 
allel stripes; 21 regul ar arrays of dropsjISE an( j irregularly 
branched structures.! 12 ! 13 ! 28 ! ^ The occurrence of these more 
complicated patterns is highly sensitive to the particular ex- 
perimental setup and the system parameters. 

Several groups employ this type of wet evaporative de- 
position as a non-lithographic technique for covering large 
areas with regular arrays of small-scale structures, such as, 
e.g., concentric gold rings with potential uses as resonators 
in advanced optical communications systems 30 or ordered 
arrays of cyanine dye complex micro-domes employed in 
photo-functional surfaces. 31 Often the patterns are robust and 
can be post-processed, e.g., to create double-mesh struc- 
tures by crossing and stacking two ladder filmsJ 2 ^ A number 
of investigations focus on deposition patterns resulting from 
more complex fluids, such as phase separating polymer mix- 
tures;^ solutions of DNA,EE3 collagen,^ liquid crystals,^ 
dye molecules,! 16 ^ 31 ! 37 ! dendrimers,^ car bon n anotubesj22E§l 
and graphene;^ and biofluids like blood.^El The latter has 
medical implications as it is thought that one may learn how 
to detect some illnesses by simple evaporation experiments on 
small samples.!^ 4 -! 

Surveying the literature, one finds that deposition patterns 
consisting of regular stripes are a rather generic phenomenon. 



They occur for many different combinations of substances^ 
and in a range of different experimental setups that allow for 
slow evaporation. To better control the contact line motion, 
various techniques are employed. One may distinguish be- 
tween passive and active experimental set-ups. In the pas- 
sive set-up, the solution or suspension is brought onto the 
substrate and left to evaporate. Examples include (i) the so- 
called "meniscus technique" (where a me niscus w ith a con- 
tac t line is created), e.g., in a sphere-on-flat) 22 ! 23 ^ or ring-on- 
flatUZH 6 geometry, (ii) the deposition of a single large drop 
onto a substrate 3 5 and (iii) the deposition of flat films onto a 
substrate using spin-coatingj 47 ! 48 ! These 'passive' set-ups are 
mainly controlled via the temperature, the partial pressure of 
the solvent, and the solute concentration. 

The active set-ups involve an additional control parameter 
that can often be better adjusted than those in the passive set- 
ups. An example is a set-up similar to blade coating where a 
controlled continuous supply of solution is provided between 
two glass plates. The upper plate slides backwards with a con- 
trolled velocity and in this way maintains a meniscus-like liq- 
uid surface where the evaporation takes place and the patterns 
are deposited.^!! The deposition patterns are found to depend 
on the plate velocity. Other examples are a receding meniscus 
between two glass plates whose receding velocity is controlled 
by an imposed pressure gradient, 26 an evaporating drop that is 
pushed over a substrate at controlled velocity,! 2 ^! or a solution 
that is spread on a substrate by a roller that moves at a defined 
speedJ 3 -^ 

Up to this point, we have only mentioned experiments in- 
volving evaporating solutions or suspensions on solid sub- 
strates. There are two systems that are closely related: On 
the one hand there exist studies of evaporating films on a fluid 
substrate, in particular, films of a nanocrystal dispersion in 
alkanes that spread and evaporate simultaneously on the free 
surface of a polar organic fluid that is immiscible with the 
alkane. 49 The defect-free liquid substrate allows for highly 
regular periodic stripe patterns that persists over a large area. 
On the other hand there are many experiments related to the 
transfer of Langmuir-Blodgett monolayers, i.e., high density 
surfactant layers, from a trough filled with water onto a solid 
plate that is withdrawn with a controlled velocity from the 
bath. Depending, e.g., on the velocity of the plate and the 
monolayer density on the trough, the transferred monolayer 
may exhibit stripe patterns parallel or perpendicular to the 
direction of withdrawal. ! 50 " 52 ! We will return to these experi- 
ments in the conclusion and discuss in which way our results 
may be related to these rather different systems. 

Despite the large variety of experiments that produce reg- 
ular line patterns from polymer solutions and colloidal sus- 
pensions, a theoretical description of their formation has been 
rather elusive. Most authors agree that the patterns result from 
a stick-slip motion of the contact line that is caused by pin- 
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ning/depinning events ffi 2 ! 30 ! 53 ! Various reduced models have 
been developed that: relate the interaction between the con- 
tact line and the deposit that is formed, in terms of a pinning 
force and derive how this force depends on and scales with the 
experimental parameters;^! 2 ^ develop evolution equations for 
the shape of an individual deposited ring; 3 study the time evo- 
lution assuming a permanently pinned contact line r 4 " 57 ! Hu 
and Larson^ analytically obtain a flow field that is combined 
with Brownian dynamics simulations to study particle deposi- 
tion. Warner et al. ^employs a thin film model similar to the 
one we present below to describe the dewetting of a film of a 
nanoparticle suspension. They study the regime where drop 
arrays form via directed convective dewetting of the solvent 
before they subsequently dry out. A thin film model that pro- 
duces rings is presented by Kaya et al.^, however, there the 
contact line is shifted 'by hand' if a certain condition is met. 
A thin film description of an evaporating drop of a suspension 
may also be coupled to a full description of the diffusion of 
vapour in the gas phased This model is used to predict the 
dependence of the mean deposit thickness on the substrate ve- 
locity, but is not used to describe the formation of deposition 
patterns. 

In another approach, the system is described using a com- 
plex non-isothermal Navier-Stokes model, i.e., with the com- 
plete set of transport equations for momentum, energy, and 
solute and vapour concentration, thereby incorporating evap- 
oration, thermal Marangoni forces and rules for contact line 
motion. 53 Bhardwaj et further incorporate Derjaguin- 
Landau-Verwey-Overbeek (DLVO) interactions between so- 
lute particles and the solid substrate in the form of effective 
forces in the advection-diffusion equation for the solute con- 
centration. Simulations show the formation of and depinning 
from a single deposit line, but for the parameter values used in 
this study, they do not observe periodic deposits.^ 

To our knowledge, there exists no efficient model that is 
able to capture the dynamics of the periodic deposition pro- 
cess, i.e., the stick-slip character of the process. By 'efficient' 
we mean a model that allows for a numerical exploration of 
the parameter space. The full transport equations contain all 
or most of the physics, however, they are tedious to use and 
tricky when it comes to incorporating details like the wetta- 
bility and contact line motion. Models that assume a pinned 
contact line are only able to describe how a deposit forms for a 
fixed drop base, even if they are fully dynamic thin film mod- 
els. For instance, the model by Okuzono et al. 55 distinguishes 
a fluid and a gel-like part of the drop and allows one to follow 
the time evolution of the drop and concentration profiles. It 
can distinguish between final deposits of basin-, crater-, and 
mound-type. The crater-type deposits might be seen as result- 
ing from the deposition of a single ring. However, models like 
this one that fix the drop base are not able to capture the de- 
position of multiple rings or of a regular line pattern in planar 



geometry. Note also that a first fully dynamic model exists for 
the related phenomenon of the transfer of Langmuir-Blodgett 
surfactant monolayers from a bath onto a solid substrate.^ 
There, the stripe formation is related to a phase transition in 
the surfactant layer that results from a substrate-mediated con- 
densation effect. We discuss similarities and differences of the 
model by Kopf et al.^ and our thin film model for a receding 
contact line of an evaporating solution in section [2] 

The aim of the present work is to analyse in detail the re- 
sults from a dynamic close-to-equilibrium thin film model that 
allows for a generic description of a receding contact line of a 
solution in the passive geometry, i.e., without an externally im- 
posed velocity. We show that the model is able to capture the 
dynamics of the stick-slip motion of the receding three-phase 
contact line of a colloidal suspension or polymer solution on 
a solid substrate. For a range of parameter values, as the sol- 
vent evaporates, the model predicts the deposition of regular 
and irregular patterns of lines that are parallel to the receding 
contact line. The model allows us to elucidate generic mecha- 
nisms that play an important role in many of the experimental 
systems that are studied. Furthermore, we are able to deter- 
mine the influence of varying the basic model parameters on 
the deposition patterns, in particular, on the spatial period, am- 
plitude, and morphology of the deposited lines. In the present 
work we restrict our attention to one-dimensional deposition 
patterns, i.e., regular and irregular line patterns. The limita- 
tions of this approach are discussed in the conclusion. 

This paper is laid out as follows: In Section|2]the thin film 
model is introduced in the form of evolution equations for 
the film thickness and concentration profile, that are obtained 
from making a long-wave approximation. The model is dis- 
cussed in the context of related models in the literature. In 
Section [3] results from time simulations are presented. It is 
shown that solely having a viscosity that diverges at a critical 
solute concentration is sufficient to trigger a self-organised pe- 
riodic pinning-depinning process that results in the deposition 
of regular line patterns in a well defined region of the parame- 
ter space spanned by our non-dimensional control parameters, 
i.e., evaporation rate, mean solute concentration, Peclet num- 
ber, parameters related to wettability, solvent vapour chemical 
potential, and the intrinsic viscosity. The final Section[4]gives 
our conclusions and an outlook on future work. Two appen- 
dices discuss details of our numerical scheme (Appendix [A) 
and the measures we use to quantify the obtained line patterns 
(Appendix|B]), respectively. 

2 The thin film model 

We consider a thin film of an evaporating partially wetting 
nanoparticle suspension (or polymer solution) on a flat solid 
substrate in contact with its vapour (see Fig. [lj. We assume 
that the densities of the solvent and solute are matched, so 
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Fig. 1 Diagram showing a liquid front that recedes due to 
evaporation and convection. The front has varying velocity v(f ) and 
the liquid film thickness profile is given by the function h{x,t). In 
the bulk film the profile c(x,t) is the vertically averaged solute 
concentration profile and outside the liquid film, the solute layer 
thickness is given by h p (x,t) = c(x,t)h(x,t). 



there is little or no sedimentation of the solute within the sol- 
vent, and slow evaporation so that we may neglect the depen- 
dence on the vertical coordinate of the solute concentration 
within the liquid film. Assuming that all surface s lopes are 
small, one may employ a long-wave approximation!*^!^ an( j 
derive two coupled evolution equations for the film thickness 
profile h(x,t) and the concentration of the solute c(x,t), which 
can be written in compact form as a pair of continuity equa- 
tions for the solution and the solute, respectively: 



d t h = -d x j c (h,c)-j e (h), 
d,(ch) = -d x [cj c (h,c)+jn(h,c)] , 



(1) 

(2) 



The equation for the solution, i.e., the film thickness evolu- 
tion equation ([TJ, contains a non-conserved evaporative source 
term and a conserved convective transport term. The equation 
for the solute i.e., the evolution equation Q for the effective 
solute layer thickness /i p — ch only has a conserved dynam- 
ics, which is made up of two terms: the first is a convective 
transport term and the second describes transport due to dif- 
fusion. These coupled convective, evaporative, and diffusive 
fluxes are given by 



-Qo(h,c)d x p(h) 



3T]( C ) 



Mh) = £\p(h)-n P ], 

r 

jd(h,c) = -Q d (h,c)d x c = -D(c)hd x c. 



d xP (h), (3) 

(4) 
(5) 



The first term on the right hand side of Eq. ([T]) is the con- 
served part of the dynamics, where j c [given in Eq. <|3j] is 
the total horizontal volume flux (integrated over film thick- 
ness). This flux is driven by the pressure gradient. The pres- 
sure p is discussed below. The convective mobility function 
Q c (h,c) — h? /3f}(c) in j c [Eq. |3j] results from the Poiseuille 
flow of the liquid in the case of a no-slip boundary condition 
at the substrate and the free upper surface. The dynamic vis- 
cosity, 17(c), exhibits a strong nonlinear dependence on the 
local solute concentration and obeys the Krieger-Dougherty 
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0.2 0.4 d) 0.6 

Fig. 2 Dependence of the viscosity on the scaled solute 
concentration <j) = c/c c , as described by the dimensionless form of 
the Krieger-Dougherty law in Eq. ( |16) . The solid line is for the 
exponent v = 1.575, whereas the dashed one is for v = 2. The linear 
Einstein relation is given as dotted line. 
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17(c) = 170 1- 



(6) 



where 770 is the dynamic viscosity of the pure solvent. Note 
that the solute bulk concentration field c(x,t) is a dimension- 
less volume fraction concentration and c c is the critical value 
at random close packing, where the viscosity diverges, i.e. 
rj (c) — > 00 when c — > c~ . For hard-spheres c c = 0.63. The vis- 
cosity concentration dependency Eq. (|6]i is illustrated in Fig. [2] 
In general, the precise value of the exponent v depends on 
the type of suspension employed. Here we mainly consider 
particles that have no net attractive forces between them and 
only have excluded volume interactions, i.e., they interact via 
a hard-core repulsion when coming into direct contact. For 
such particles, values for v between 1.4 and 3 are discussed 
in the literature, depending on their shape. 6 ^ The exponent v 
is sometimes written as v = [rj)c c , where [rj] is the intrinsic 
viscosity, see,Larson 66 defined by 



tj = hm — . 

c^Q TJ C 



(7) 



For spherical particles [tj] = 2.5, r esulting in v = 1 .575. Other 
thin film models use v = 2J 59 I 68 I For solute particles with net 
attractive forces between them, values for v as low as 0.13 are 
reported. 69 Depending on the particular system, the transition 
at c c is either referred to as jamming or gelation.'^^ 9 -' Here 
we fix v = 1.575, corresponding to particles that only interact 
via excluded volume interactions, in order to gain an under- 
standing of these types of systems and also to determine the 
regions of qualitatively different behaviour in the phase plane 
spanned by the evaporation rate and concentration, defined be- 



low. Additional investigations show (see Section 3.7 1 that the 
formation of periodic deposits is somewhat more pronounced 
for smaller v, which corresponds to more strongly attracting 
solute particles. 
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Fig. 3 (Color online) The non-dimensional disjoining pressure 



n(h), given in Eq. 1 15 I, for various values of %, as given in the 



legend. The horizontal straight lines indicate various values of the 
non-dimensional chemical potential M, given in Eq. l |19} . The blue 
solid lines for the disjoining pressure and chemical potential 
correspond to our 'standard' values of % = 1.085 and M = —0.003, 
respectively (see the end of Section[2](. The other lines correspond to 
parameter values employed in Section ing] below. The blue stars 
denote the stable film heights: hi < hi for the standard values. 



The convective flow is driven by the gradient of the pressure 

p(h) = - r d xx h-n(h), (8) 

where the first term is the Laplace pressure (y is the surface 
tension) and the second term is the disjoining (or Derjaguin) 
pressure 

S F do — h 



n(h) 



2« 2 



- exp - 



(9) 



h 3 Iq l 
that models wettability effects for a partially wetting 
fluid ;IS5|2QED see Fig. [3] for a graphical representation of its 
non-dimensional form. Iq is the Debye length, do is a molec- 
ular interaction length, 5 LW = — AjYlnd^ and S p < are the 
apolar and polar spreading coefficient, respectively, and A < 
is the Hamaker constant. We expect qualitatively similar be- 
haviour for other combinations of stabilising and destabilizing 
terms in il(A).™] 

The second term on the r.h.s. of Eq ([TJ represents the 
non-conserved part of the dynamics and models evaporation, 
where j e [given in Eq. Q] is the evaporative flux density at 
the free surface of the film. Here, we assume that the system 
is close to equilibrium and that the vapour is near to saturation 
and so evaporation is slow. In this limit evaporation with a rate 
j3 is driven by the difference of the scaled pressure p /p and 
the chemical potential of the ambient vapour jU. E3ZH Latent 
heat effects may be neglected, and the density p is assumed to 
be equal for the solute particles and the solvent. 

The first and second terms on the right hand side of Eq. Q 
model convective and diffusive transport of the solute, respec- 
tively. Since the solute is passively advected by the fluid, the 



convective flux is given by cj c . For the diffusive flux [given 
in Eq. ([5]l] we assume Fick's law and we set the diffusive mo- 
bility to be Qd{h,c) — D(c)h, where D(c) is the concentration 
dependent diffusion coefficient, and we employ the Einstein- 
Stokes relation 



D(c) 



6nror]{cY 



(10) 



where kg is the Boltzmann constant, T the temperature, and ro 
the solute particle radius. 

To facilitate a comparison of our results with the literature 
we employ a scaling identical to that employed by Lyush- 
nin et alM^ and,Frastia et alM* namely: time, horizontal x- 
coordinate, and film thickness scales are 



x = 



and 



respectively, where 
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(12) 



Using this scaling, Eqs. ([T]) and Q are brought into the fol- 
lowing non-dimensional form 



d t h = d x 
d,((j)h) = d x 

where 



h 5 



W) 



d x p(h) 
d x p(h) 



■Q \p(h)-M 
h 



PeT}(0) 



d x <j) 



(13) 
(14) 



p(h) = -d xx h-n(h) = -d^h-hT 5 + exp(-xh). (15) 



Note that starting with Eqs. ( 13 1 and ( 14 1 in the remainder of 
the paper, the symbols h, h p , t, x, p, 7], and 11(h) stand for the 
non-dimensional quantities whereas up to this point they de- 
noted the dimensional quantities. However, for the scaled con- 
centration we introduce <p =c/c c . The diffusion coefficient is 
expressed as [Perj(0)] _1 , where Pe is the Peclet number and 
the dimensionless viscosity function is 



77(0) = (1 -0)- V - 



(16) 
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The dimensionless numbers in Eqs. ( 13 1-( 15 1 are defined as 

(17) 
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(18) 

(19) 
(20) 



where the evaporation number, Qq, is the ratio of the time 
scales of convection and evaporation of a film without solute 
and the reciprocal Peclet number, Pe~ 1 , is the ratio of the time 
scales of convection and diffusion. 

There exist a number of models in the literature that are 
similar to that defined in Eqs. ( 13 1 and ( [T4) i. One is used to 
study macroscopic particle-laden film flow down an incline.^ 
There, the solvent is non-volatile and changes in the particle 
concentration result from the settling of particles due to grav- 
ity. As the particles in the model of Cook et al. 68 are assumed 
to be large, no diffusion is included in the model, nor is wet- 
tability; the advancing of the contact line is facilitated by a 
precursor film of imposed height. Another example is the case 
of a surface-passive solute studied by Warner et al.^ [their 
Eqs. (54) and (55)]. They use a different Derjaguin pressure 
that, however, also models partially wetting liquids. The main 
difference is in the evaporation model that they use, which in- 
corporates a vapour recoil effect. They also model particle 
diffusion in a manner that is independent of solute concentra- 
tion. Both, Cook et al. 68 and Warner et alJ^ use a Krieger- 
Dougherty law to model the dependence of the viscosity on 
the concentration. This is also done by Craster et alM^ where 
the spreading and retraction of evaporating droplets containing 
nanoparticles is studied with a thin film model that involves a 
structural di sjoining pressure. There are some similarities to 
other modelsl 54 ! 55 ! 57 ! where, however, the contact lines are kept 
pinned. Since moving contact lines are an intrinsic part of the 
solute deposition process that we study here, we do not refer 
any more to these other approaches. 

The evaporation model that we use is valid close to equilib- 
rium where the film may be considered isothermal. The evap- 
oration is limited by the kinetics of the phase transition (or by 
the boundary layer transfer, but not by the diffusion of sol- 
vent vapour in the gas phase) in the contact line region, which 
is influenced by the effective molecular interactions, i.e., the 
saturated vapour pressure depends on the disjoining pressure 
and curvature in the ma nner e mployed in studies of evapo- 
rating films and drops J 74 " 76 ! 78 ! 79 ! Our evolution equation q 1 3 p 



Note, that one may also obtain our evaporation model by tak- 
ing the isothermal limit of the models by.Ajaev^, Rednikov 
and Colinet 81 For a further discussion see.Todorova et al.^i 
In the present work, we restrict our attention to line pat- 
terns deposited through evaporative dewetting, i.e., we as- 
sume a one-dimensional geometry, as sketched in Fig. [T] It 
is known from the experiments that line patterns are not al- 
ways transversally stable. We discuss this point in relation to 
our results in the conclusion. We start our numerical com- 
putations with an initial condition that corresponds to a spa- 
tially one-dimensional semi-infinitely extended film of con- 
stant thickness that is connected to a thin precursor layer. The 
thickness of the precursor film corresponds to the equilibrium 
height where the Derjaguin pressure and evaporation balance. 
The liquid in the film is a suspension with a constant solute 
concentration. We discretize the non-dimensional equations 
13} and([l4} over a finite domain x £ [0,L] where the dewet- 



reduces to the model by,Lyushnin et alM^in the limit <j)Q — > 



ting front is located close to the boundary of the domain at 
x = 0, see Fig. [T] The details of the implementation including 
the procedure we use for shifting the spatial computational do- 
main at certain times are given in Appendix [A] 

Initially, we employ a set of values of our dimension- 
less parameters that we will refer to as the standard config- 
uration: We set the dimensionless chemical potential, M = 
—0.003, and the dimensionless disjoining pressure parameter 
X = 1.085, which are the same values as used by.Lyushnin 
et al.^ The other parameters are modified to accommodate 
the effect of the solute. In particular, for the initial results pre- 
sented here, we use the evaporation rate — 4.64 x 10~ 7 ; 
the scaled constant bulk concentration, <po = 0.41, that enters 
the simulation as a given initial value; the reciprocal Peclet 
number, Pe~' = 0.0003; and the exponent of the Krieger- 
Dougherty law, v = 1.575, as discussed above after Eq. (fT). 

The main control parameters in our simulations are Qq and 
00- These two parameters span the plane within which we de- 
termine regions where periodic patterns and other characteris- 
tic deposits occur. Then, for selected fixed parameter combi- 
nations (Qo,^o), we vary one of the parameters Pe -1 , M, %, 
and V while the remaining parameters are kept fixed, in order 
to determine how sensitive the behaviour of the system is to 
variations in the value of these. 



3 Results 

3.1 Dynamics of the evaporative dewetting front 

Our investigations are mainly based on time simulations of 
the front motion and the deposition process. We start with a 
semi-infinitely extended liquid film that coexists with an ul- 
trathin precursor film. The initial dewetting front has a step- 
wise profile, but is smoothed by capillarity during the very 
first time steps. The concentration is initially set to a uniform 
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Fig. 4 Sequences of snapshots of film thickness profile h(x, t) (top 
panel), mean concentration <p(x,t) (middle), and evaporation flux 
je(x,t) (bottom), which illustrate the pinning-depinning dynamics. 
The individual profiles in each sequence are plotted starting at the 
same time, with a time-increment At = 10 9 between each successive 
curve, and are displayed with periodically repeated line styles (solid, 
dashed, dash-dotted, dotted, solid, ...). All the model parameters are 
set to the 'standard values', which are defined at the end of Sec.|2| 



value (j)(x,0) = 0rj. The chosen value of the vapour chemical 
potential ensures that the front recedes by evaporation and/or 
convection. As the front recedes, it deposits part of the solute 
in an initially smooth layer. Then, in most cases, the front set- 
tles after some transient into a different type of regular motion. 

The time-dependent behaviour of the system is well known 
in the case without solute (0o = 0).^ There, the front profile 
always converges to a constant shape that moves with constant 
velocity, i.e., the front motion is stationary. In this situation, 
one may still distinguish between two qualitatively different 
limiting cases: (i) convection-dominated and (ii) evaporation- 
dominated dewetting. They are found for small and large val- 
ues of the evaporation number i2o, respectively. In case (i) 
the convective motion maintains a capillary ridge despite the 
ongoing evaporation whereas in case (ii) convection is much 
slower than evaporation and there is no capillary ridge. 

In the presence of a solute, the situation is more complex 
and there exist extended regions in parameter space where 
no stationary front motion is found. Instead, the receding 
front shows an unsteady motion with periodically changing 
front velocity and shape. An example of such a dynam- 
ics is shown in Figs. |4] [5] and [6] These figures illustrate 
various aspects of the periodic pinning-depinning dynamics 
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Fig. 5 Snapshots from the deposition process at two characteristic 
stages: evaporation-dominant (top) and convection-dominant 
(bottom). The system parameters are as in Fig. [4] The profiles in the 
top and bottom panel correspond to the 4th and the 1 1th profile in 
time in Fig. [4] respectively. Note the capillary ridge behind the front 
in the lower plot, which is absent in the upper 
evaporation-dominated stage of the cycle. 

of the front that is related to a periodic transition between 
evaporation- and convection-dominated regimes of the front 
motion. Fig. |4] shows snapshots of the film height h(x,tj), 
concentration <p(x,tj) and evaporation flux y e (jc,f,) profiles. 
They are equidistant in time and cover just over one period 
in time. The pinned (evaporation-dominated) regime is char- 
acterised by densely spaced profiles when the front is at about 
x = 4.275. During this stage the front moves very slowly, the 
height profile has no capillary ridge and the evaporative flux is 
rather localised. In contrast, during the depinned (convection- 
dominated) stage of the cycle, the profiles are sparsely dis- 
tributed (i.e., the front is fast), the convective motion supplies 
enough solution to maintain a capillary ridge and therefore the 
evaporative flux is spread over a wider x -region. Note that the 
maximal j e (x,t) is actually larger during the depinned regime. 

Transitions between the two regimes can be explained as 
follows: In the early stage of the convection-dominant phase 
the front moves relatively fast, although the evaporation flux, 
j e (x,t), is still greatest in the contact line region. This re- 
sults in an increase of the local concentration <j)(x,t) and, in 
consequence, in a strongly nonlinear increase of the viscos- 
ity in this region that suppresses the convective motion of the 
front. As a result, the front slows down and if the local so- 
lute concentration reaches random close packing (0 (x. t ) — > 1) 
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Fig. 6 (Color online) Space-time plots of the nanoparticle deposit 
thickness h p (x,t) (top) for the standard values of the system 
parameters, where the deposited lines have characteristic heavy right 
tail (see Fig. 1 1 1 and a contour plot of the film height h(x,t) 



(bottom) in the rectangular region marked in the upper plot. There 
are 24 equally spaced level contours in the interval of heights 
[1.2,3.25] which allows the steep dewetting front and capillary 
ridges to be resolved. 





IS 



Fig. 7 (Color online) Space-time plots of the nanoparticle film 
height h p (x,t) for 2 characteristic cases: the upper panel shows 
deposit lines with a dominant left 'tail' (Qq = 4.64 x 10~ 6 , other 
parameters are standard). The regular line pattern corresponding to 
this situation is shown in Fig.[TT| The lower panel shows a 
long-period large-amplitude pattern (<fo = 0.31, other parameters are 
standard). A regular line pattern for the close, but slightly different 
value 0o = 0.3016 is shown in Fig.jTTJ 



the convective motion in the contact region stops completely. 
Thus, the suspension becomes locally jammed. In the case 
of a polymer solution the transition may be referred to as a 
local gelling transition (as employed, e.g., in the piece-wise 
model by.Okuzono et alP^) With the convective motion ar- 
rested, the front only moves by evaporation. As the resulting 
front velocity can be orders of magnitude slower than during 
the convective motion; the front seems to be pinned. This is 
clearly visible in the space-time plots of Figs. [6] and [7] The 
typical front shape at this stage is then monotonic, and no cap- 
illary ridge exists [see the solid line in Fig. |5ja)]. During this 
phase of slow evaporative motion, the front effectively leaves 
deposits of the highly concentrated solute. As a result, the 
concentration and therefore viscosity in the contact line re- 
gion decrease, allowing the convective motion to start again. 
The front speed is then much larger than in the evaporation- 
dominated phase, and the front seems to depin. The typical 
front shape at this stage is non-monotonic, due to the presence 



of a capillary ridge [see the solid line in Fig.[5|b)]. Note that 
the width of the region at the front where the concentration 
of the solute is locally increased varies over time. The ex- 
tent (along the x-axis) of this region is visible in Fig. [7] as the 
bright (white) region of the field hp(x, t) at the receding dewet- 
ting front. We see that in the convective regime the width in- 
creases as the solute gets concentrated and in the evaporative 
(pinned) regime the width decreases as the solute is deposited 
and left behind the front. 

The resulting periodic pinning-depinning cycle - that is per- 
ceived as a stick-slip motion - can be best appreciated in the 
space-time plot displayed in Fig. |6jb) which shows the film 
thickness profile. The contours of h(x,t) in the steep dewet- 
ting front region are so closely bunched that they appear to be a 
single thick line whose slope, dx/dt, corresponds to the veloc- 
ity of the dewetting front. It is clearly visible that the veloc- 
ity periodically changes between two rather different values. 
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Fig. 8 (Color online) A selection of typical deposit profiles, 
including the spatial region where the transient to long-time 
behaviour occurs. The panels show deposition of (a) no lines (but in 
the transient stage h p may change nonmonotonically); (b) a single 
line followed by a flat layer; (c) transient lines (whose amplitude 
decays first fast then slow) followed by a flat layer; (d) transient 
lines (whose amplitude decays first slow then fast) followed by a flat 
layer (decreased v = 1.25); (e) an intermittent line pattern [which is 
a magnification of the pattern in panel (g)]; (f) transient double lines 
(converging to regular lines); (g) an intermittent line pattern; (h) 
transient lines followed by a regular line pattern; (i) transient lines 
converging to a regular pattern of double lines (decreased 
X = 1.065); (j) a long-period pattern switching between a flat layer 
and a single line with a leading depression (decreased v = 1.25). 
The corresponding parameters for (a-j) are: Qq = (4.64, 14.68, 
4.64, 0.167, 0.147, 0.464, 0.147, 0.464, 0.7, 0.464) x 10~ 6 , 
<t> Q = (0.3343, 0.41, 0.3588, 0.41, 0.41, 0.2983, 0.41, 0.41, 0.41, 
0.498). The remaining parameters are the standard values, defined in 
section|2] The letters (a-d, g, h) are used to mark the corresponding 
regions in the phase diagram displayed in Fig.|9] 



Further contour lines to the right of the moving front show 
the periodic appearance of a capillary ridge in the phases of 
convective motion. The pinning-depinning cycle can also be 
seen in the space-time plots of the effective nanoparticle layer 
thickness h p (x, t) presented in Fig.[6ja) and|7] There, however, 
the capillary ridges are not as clearly visible because h p (x,t) 
decreases faster when moving away from the dewetting front 
into the wet region. 

The significant changes in the relative importance of con- 
vective and evaporative fluxes can be conveniently described 
by a local evaporation number £2\ oc = (rj (0) /T7o)X2o- When 
Q\ oc diverges, the front seems pinned, but actually still moves 
extremely slowly by evaporation alone, and deposits a line or 



a thick layer of solute. When Q\ oc sufficiently decreases, the 
front depins and convective motion resumes. 

We use the simulation set up described above to investi- 
gate the transient and long-time deposition behaviour. It is 
found that after the initial transient (that may involve the de- 
position of some irregular lines) various scenarios are possi- 
ble, depending on the values of our model control parameters. 
Most importantly, there exists an extended region in parame- 
ter space, where after a transient, very regular line patterns are 
deposited (see Figs. [5] [6] and [7). In these cases the pinning- 
depinning process is repeated in a regular manner. Particular 
features of the deposited pattern depend on the particular val- 
ues of the model parameters and are discussed below. The oc- 
currence of periodic deposits for an extended parameter range 
indicates that this phenomenon is robust and explains why the 
deposition of line patterns is found in a wide variety of experi- 
mental settings with various suspensions and solutions, where 
it is often described as resultin g from a regular stick-slip mo- 
tion of the contact line. | 22 | 26 | 3 ° l 

Before we discuss the various types of deposition patterns 
that we observe as the system parameters are varied, we return 
to the space-time plots in Figs. [6] and [7] that already indicate 
how the ratio of typical time scales for the convection- and 
evaporation-dominated front motion changes with increasing 
evaporation number and mean concentration. Fig. [6] shows 
(as do Figs. |4] and [5]) our standard case (described at the end 
of Sec. |2j. There, the ratio of the velocities of the dewetting 
front in the depinned and pinned regime is 45.07. Increasing 
the evaporation rate coefficient to £2q — 4.64 x 10~ 6 (keeping 
the other parameters at the standard values), the evaporative 
flux becomes stronger and the convective motion is not fast 
enough to create a large capillary ridge. As a result, the differ- 
ence between the depinned and pinned regimes is smaller and 
so the ratio of the velocities decreases to 8.16 [see Fig.|7|a)]. 
Further increasing of Qq eventually results in the deposition 
of a layer of constant thickness. If instead one decreases the 
concentration to 0o = 0.31 (keeping the other parameters at 
standard values), see Fig. [7jb), the viscosity in the convective 
regime decreases and the capillary ridges become larger. Be- 
cause there is less solute present, it takes longer to build up a 
high concentration of solute in the contact line region and to 
pin the front. This results in an increased length of the con- 
vective phase of the cycle (compare Figs.|7jb) and|6| and to an 
increase of the ratio of velocities to 80.36. A further decrease 
in the concentration prolongs the convective phase even more 
until eventually the length of the convective phase diverges (at 
finite concentration value, see below). 

3.2 Types of deposition patterns 

An extensive parameter scan in the space spanned by our con- 
trol parameters - £2q, 0q> Pe~ ! , %, M, and v - reveals a zoo of 
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Fig. 9 (Color online) Morphological phase diagram of deposition 
patterns, in the plane spanned by the evaporation number £2q and the 
bulk concentration 0o- Letters (a)-(d), (g) and (h) in the differently 
shaded areas (coloured online) indicate the type of pattern that is 
found there using the same letters as in the caption of Fig. [5] In 
particular, in the central grey area (green online) regular line 
patterns are found after some transient (simulations denoted by o), 
while outside this area (simulations denoted by *) a layer of constant 
height is deposited after a variety of transients indicated by 
subregions of different shading (colour online). As an exception, an 
intermittent pattern is found in the small (red online) region (g). 
Results obtained along the dashed vertical and horizontal lines are 
presented in subsequent figures and discussed in detail in the main 
text. The standard configuration defined at the end of Section[2]and 
shown in Figs.|4]to|6] is marked by the small filled square where the 
left vertical and upper horizontal dashed lines cross. Fig.[7]gives 
results for the two other crossing points that are marked by filled 
triangles. 

various different deposition patterns and allows one to study 
the dependence of the pattern morphology on the control pa- 
rameters. An overview of typical deposit profiles, (x), is dis- 
played in Fig. [8] The letters (a)-(j) that indicate the individual 
panels are used in the following list that describes important 
properties of the individual patterns. They are also used to in- 
dicate the corresponding regions in the morphological phase 
diagram in Fig. [9] 

(a) After a short transient during which the height of the de- 
posit can vary nonmonotonically (there can be a very 
small bump but the solution is not in jammed state), 
a layer of constant thickness is left behind the moving 
front. This is normally observed for very dilute solu- 
tions at any evaporation number. The solute concentra- 
tion in the resulting film is everywhere below the jam- 
ming threshold, i.e. the solute may still diffuse within the 
precursor film. This is in contrast to all the other cases 



where at least parts of the deposit (normally the lines) 
are in the jammed state. 

(b) A single line is deposited as the final part of a long initial 
transient before a flat layer is deposited. This represents 
the case closest to the original coffee stain effect where a 
single ring is deposited from an evaporating droplet. This 
behaviour is observed for a wide parameter range outside 
the region of periodic patterns [see (h) below] and for 
denser solutions than in case (a). Apart from in the region 
of the initial line deposit, the solute in the flat layer is not 
necessarily jammed. 

(d) A finite number of lines is deposited as part of a long 
initial transient before a flat layer of solute is deposited. 
This represents the experimental case where multiple ir- 
regular rings are deposited from an evaporating droplet. 
This behaviour occurs in bands around the region of peri- 
odic patterns. Panels (c) and (d) in Fig. |8]present deposit 
profiles from the regions with higher and lower than 
in the periodic line region (cf. Fig. [9]), respectively. 

(f) A very long transient of a periodic pattern of double lines. 
The pattern slowly evolves towards one with regularly 
spaced lines by slowly changing the distance between 
the two lines in the respective pairs. The initial tran- 
sient deposit is a flat layer and the size and shape of the 
first line is not significantly different from the following 
lines. This pattern is observed for very low solute con- 
centrations and for somewhat small evaporation rates in 
the region of periodic lines close to its boundary. 

(g) An intermittent pattern (a magnification of a portion of 
this is displayed in Fig. [Hfe)). A rather irregular line 
pattern is deposited in an intermittent manner, i.e., there 
are times when a nearly flat layer is left behind the 
front which alternate in a non-periodic way with episodes 
where line patterns are deposited. These line patterns are 
not periodic but exhibit a typical timescale; see the mag- 
nification in Fig.|8je). This behaviour is found in a very 
narrow band between the region of periodic patterns and 
the band where single or multiple lines form (cf. Fig.|9]l. 

(h) After a short monotonous transient, the deposition of 
lines starts with a large first line followed by a small num- 
ber of transient lines. The transient lines in Fig.[8jh) are 
of increasing period, however, the case of decreasing pe- 
riod is also observed. The line deposition rapidly con- 
verges to a regular periodic state. For instance, Fig. [HJh) 
shows a case where a rather high initial line terminates 
in a shoulder that is then followed by lower lines that 
slowly increase in height and converge to the truly peri- 
odic line pattern. This profile is at the standard values of 
our control parameters; see also the results in Figs.[5]and 
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[6] This behaviour is found in a large region of the space 
spanned by our main control parameters (Fig. 15}. The 
region may shrink and disappear if other parameters are 
strongly changed away from their standard values, such 
as, e.g. increasing the diffusion coefficient leads to the 



the parameters of our numerical solvers; see Appendix A. 3 



pattern formation being suppressed (see section 3.5 1 



(i) Stable double lines can be observed for values of the pa- 
rameter x in the disjoining pressure that are less than our 
standard value. For the case displayed in Fig. [8|i), we 
observe that a large initial line is deposited followed by 
some irregular lines with monotonically increasing line 
distance. However, unlike the transient double lines in 
case (f), here the distances adjust till a regular pattern of 
double lines is deposited in a stable manner. 

(j) After the initial transient deposition of a finite number of 
lines it at first seems that a region of flat film follows. 
However, the fiat film actually becomes a long-period 
pattern that switches between a fiat layer and single lines 
with a leading depression. This pattern was observed for 
a lower value of the viscosity exponent, v = 1.25, that 
corresponds, e.g., to a suspension of non-spherical parti- 
cles. A similar pattern was observed at the boundary of 
region (h), where 0o = 0.31 and £2q is small. In this case, 
switching between a flat layer and multiple lines occur. 

Note, that in some cases where the concentration is small, the 
concentration in the 'valleys' between the periodic lines is suf- 
ficiently low that the solution is not jammed in these regions; 
an example is given in Fig. [8jf). 

As our parameter space is 6 dimensional and the equations 
of the model (13i and (14i are highly nonlinear, our list of 



typical patterns must certainly be incomplete. As our inves- 
tigation is numerical in nature, we are unable to determine 
to which state some of the very long initial transients con- 
verge. Below, we discuss the various transitions between the 
observed patterns on changing a control parameter. This al- 
lows us to speculate with more confidence which further types 
of patterns might be expected. 

Although one is easily able to qualitatively classify the 
transient patterns, a quantitative analysis is cumbersome and 
therefore is not pursued here. Discarding the initial transients, 
we can nevertheless distinguish several types of deposits: (a)- 
(d) flat layers; (f), (h) regular periodic lines; (g) intermittent 
line patterns; (i) periodic arrays of double or multiple lines; 
(j) periodic switching between a depression-line combination 
and a flat layer. Our main aim here is to analyse the periodic 
line patterns that are observed in region (h) of Fig. [9] This is 
done in the following sections, where we analyse the depen- 
dence of the line morphology on selected control parameters, 
while keeping the other parameters fixed. We have carefully 
checked that the patterns are robust by using various values of 




Fig. 10 (Color online) Measures characterising the regular lines as a 
function of tj>Q, for fixed standard Qq = 4.64 x 10 -7 (blue line with 
*) and for fixed £2 Q = 4.64 x 10~ 6 (green line with o).EllThe other 
parameters are fixed at their standard values. Selected corresponding 
line patterns are displayed in in Fig. 
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3.3 Regular patterns in the (Qq, 0o) plane 

Having discussed the dynamics of the deposition process (sec- 
tion |3. 1| > and the main types of transient and long-time deposi- 
tion patterns (section [3~2| i, we now embark on a more detailed 
analysis of the morphology of the regular line patterns and its 
dependence on the location in the parameter plane (i2o,0o)- 
On the basis of our exploration of this parameter plane we de- 
tected regions of characteristic deposition patterns - see Fig.[9j 
The sub-regions marked by letters in Fig. [pjrelate to the pat- 
terns under the same letters in Fig. [8] and the list of typical 
patterns in Sec. 3.2 The following analysis is based on a large 
number of long-time simulations (for details see Appendix [A) 
in the region (h) of periodic lines, cf. Fig. |9j Excluding the 
initial transient, we take a sequence of TV regular deposition 
periods (lines), where 10 ;$ TV ;$ 100, depending on the spatial 
period of the deposit and the required CPU time. We process 
the profiles and extract measures that characterise the individ- 
ual lines and their arrangement. In particular, we obtain the 
amplitude a, spatial period P, relative width K = 2a /P (where 
a is the standard deviation), skewness 73, kurtosis 74, excess 
cross-sectional area of the lines A, the mean deposit height 
hp, and the minimal height of the deposit between the lines 
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Fig. 11 (Color online) Morphology changes of the line pattern as 0o 
is varied, corresponding to Fig.[l0] Left column: For fixed standard 
£2o = 4.64 x 10~ 7 and, from the top to the bottom, 
00 = 0.3016, 0.34, 0.41 , 0.475, 0.494, 0.5, and other parameters are 
standard. Right column: For fixed Qq = 4.64 x 10~ 6 and, from the 
top to the bottom, fa = 0.389,0.4,0.41,0.46,0.49,0.499, and other 
parameters are standard. 



hp. Definitions of all these measures are given in Appendix|B^ 
Our investigation shows that these quantities strongly depend 
on both the evaporation number Qo and the concentration (j>o- 
We focus on two vertical and two horizontal cuts through the 
parameter plane (Qo, <po) that are indicated by dashed straight 
lines in Fig. [9] 

First, we vary the concentration 0o for two evaporation 
numbers, Q = 4.64 x 1CT 7 and Q Q = 4.64 x 1CT 6 , while the 
remaining system parameters are fixed at our standard values 
(see end of section^. Fig. 10 presents the line characteristics 
for both cases and Fig. [TT] shows a number of corresponding 
deposit profiles. 

Most of the line characteristics behave qualitatively similar 
for the two cases. We first focus on Qo = 4.64 x 10~ 7 and 
then point out the differences. On increasing <po from a re- 
gion where there is no deposition of periodic lines, one first 
encounters large amplitude almost solitary peaks that are sep- 
arated by very large distances. On further increasing Qo, the 
amplitude first hardly changes and later decreases. The period 
rapidly decreases while the relative line width and skewness 
increase almost linearly with a slight drop in skewness at very 
high concentrations <po ~ 0.49. For higher Qo, the deposit pat- 
tern almost turns into a constant thickness layer, with a small 
amplitude harmonic modulation. Finally, the amplitude goes 
to zero (at finite period) at the upper border of region (h). For 
Qo = 4.64 x 10~ 6 , the evaporation is stronger and one must 
go to a higher Qo than when Qq = 4.64 x 10~ 7 , to see periodic 
deposition of lines. The period of the lines does not seem to 



diverge at this border and the amplitude of the lines is much 
smaller. Actually, the amplitude first increases with increas- 
ing Qo and has a maximum well inside the region of regular 
lines before it decreases towards zero at the other border of the 
region. Interestingly, in contrast to the former case, the skew- 
ness changes from negative to positive values, i.e., the individ- 
ual lines change their morphology from having a 'tail' away 
from the receding front to having a tail in the direction towards 
the front, as is always the case for Qq — 4.64 x 10~ 7 . In both 
cases the excess area greatly decreases, the minimum height 
of the deposited valleys between lines greatly increases, and 
the period and amplitude decrease as the boundary of region 
(h) at higher Qo is approached. The average height of the de- 
posit linearly increases with the concentration as expected but, 
interestingly, the two linear dependencies for the two evapora- 
tion rates are shifted with respect to each other. This aspect is 



further discussed in Appendix A. 2 




Fig. 12 (Color online) Measures characterising the regular lines as a 
function of £2q, for fixed standard Qq = 0.41 (blue line with *) and 
for fixed Qq = 0-31 (green line with o)JS3The rest of the parameters 
are fixed at the standard values. Selected corresponding line patterns 



are displayed in Fig. 13 



Second, we vary the evaporation number Qq for two con- 
centrations <po = 0.41 and (j>o = 0.31, while the remaning pa- 
rameters are fixed at our standard values (see end of section[2]). 
Fig. [12] presents the line characteristics for both cases and 
Fig. 1 1 3 1 displays a number of the corresponding deposit pro- 
files. 

We first focus on the case Qo = 0.41 and then point out the 
differences found for the lower concentration value. Increas- 
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Fig. 13 (Color online) Morphology changes of the line pattern as 
Qq is varied, corresponding to the results in Fig. [12] Left column: 
Fixed standard <^)o = 0.41 and, from the top to the bottom, 
A) = (0.1598,0.1667,0.3162, 1,5.274,5.663) x lO" 6 . Right 
column: Fixed 0q = 0.31 and, from the top to the bottom, 
A) = (1. 1.137, 1.468,2.15,3.162,4.083) x 10~ 7 . 



ing Qq, one moves from the narrow region (d) of transient 
multiple lines (see Fig.[9]i, passes through a very narrow band 
of intermittent line patterns (g) followed by the region (h) in 
which we observe the regular line patterns. For the lowest val- 
ues of Qo m (h)> the patterns have a relatively small period 
and a small but non-zero amplitude. The strongly anharmonic 
peaks are skewed to the right with their tail pointing towards 
the wet side. On increasing Qq the period increases. The am- 
plitude, however, first increases and then decreases, until at a 
certain threshold, the pattern ceases to be periodic and we ar- 
rive in the narrow border region (c), where only a finite num- 
ber of lines are deposited. Correspondingly, the relative width 
of the lines decreases as one moves deep into the region of 
periodic lines, where the lines become more peaked. Remark- 
ably, the skewness changes sign here, i.e. as Qq increases the 
tail of the lines shifts from pointing towards the receding wet 
film, to pointing away. 

The second cut with varied Qq is at fixed 0o = 0.31 and the 
standard values for the remaining parameters. For small Qq, 
the onset is similar as for the previous cut at 0o = 0.41. How- 
ever, on increasing Qq, one finds large-amplitude long-period 
line patterns. The period seems to diverge in a manner similar 
to that found when increasing 0o for fixed Qq = 4.64 x 10~ 7 



(see Fig. 10 1. Note that in Fig. 12 the average deposit thick 



ness changes as Qq is varied, even though the initial concen- 
tration 0o is fixed. The deposited average thickness for more 
dilute solute concentrations can be higher than when a denser 
solution is used, as long as the evaporation rate is higher. This 



effect stems from the influence the moving evaporation front 
has on the concentration field in the thick liquid film. This is 
related to our boundary conditions and is further discussed at 



the end of Appendix A. 2 



Note, that a change in the sign of the skewness as described 
above, was also observed in experiments on nanoparticle sus- 
pensions 83 and can be explained as follows: For smaller val- 
ues of Qq and/or higher 0o> the capillary ridge is large and ac- 
commodates a large amount of solute as it recedes. When the 
front pins, the capillary ridge evaporates and the solute con- 
tained in the ridge is deposited in a thick tail pointing towards 
the receding front. When the liquid front depins, it carries 
much of the solute within the liquid front away with it, which 
results in a further drop in the deposition thickness (seen as 
the final shoulder of the tail). For higher values of Qq the 
capillary ridge is smaller and so the tail pointing towards the 
receding front is smaller. On the other hand a negative skew- 
ness is typical for higher Qq and (possibly) low 0o- There, one 
encounters smaller tails towards the receding front and also a 
pronounced tail away from the front. In this situation, it takes 
some time (during which the front travels some distance) for 
the solute concentration to build up at the front. During this 
time the front gradually slows down and deposits the solute 
with growing thickness until local jamming occurs. Then the 
front pins and starts to deposit a line. The resulting asymmetry 
is seen as negative skewness. 

3.4 Onset of formation of periodic deposits 

To further clarify how the onset of pattern formation occurs we 
next focus on the small band in the parameter plane (Qq, 0o) 
that bounds the region of periodic deposits. From inspecting 
Fig.|9]and the line deposition descriptions outlined in the pre- 
vious section, one clearly sees that the onset of the formation 
of periodic deposits may occur through a number of different 
transitions. 

The most intricate transition is the one that involves the in- 
termittent deposits shown in Fig. |HJg) and (e), and occurs in 
the narrow region (g) of Fig. [9] We find that this behaviour is 
very sensitive to computational details (which is not the case 
for the regular line patterns), a fact which bolsters our opinion 
that the intermittent line patterns represent a 'chaotic depo- 
sition' . Taken in this context, the occurrence of periodic de- 
posits of double lines (Fig.|8jf),(i) would be related to the peri- 
odic deposition of single lines through a period doubling bifur- 
cation that occurs when changing the relevant control param- 
eter. Both effects point towards the presence of chaos in the 
system, as they are elements of the intermittency and period- 
doubling route to chaos, respectivelyiS] However, to investi- 
gate these effects further, simpler asymptotic models for the 
moving material-depositing front need to be constructed that 
are numerically less challenging. 
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0.0197 


0.0299 


0.0596 


0.0102 



Table 1 The measured hysteresis intervals, A. d ^ s and 4 re i, are listed for the various transitions between the long-time deposition of periodic 
patterns and flat layers that occur along the straight lines in Fig. [9] In these intervals different deposition patterns may be obtained, depending 
on initial conditions. The first row indicates which parameter is held fixed whilst the second row states which parameter is varied and which 
border of the region of periodic deposits is being considered. The third row gives the absolute hysteresis interval A. d b s in terms of the 
parameter that is being varied, whereas the fourth row gives the relative interval size 4 re i = A a \, s /A, where A is the corresponding size of the 
entire region of periodic deposits in the (logj () .Qo,0o) pl ane = 1-8239 in log 1() f2o coordinate and A = 0.225 in 0g coordinate). Note that 
4 a b s is measured in terms of the particular parameter that is being varied (either log 10 £2q or 0o)> i- e -> it represents a length in the 
semilogarithmic phase plane Fig.[9] Note, that we do not include results for fixed £2q = 4.64x 10~ 7 at low There the line period becomes 
very large and long pieces of nearly flat layers occur between the lines. This results in rather expensive computations. However, our limited 
results indicate that there is no hysteresis in this case. 



In general, it has proved to be difficult to detect the exact lo- 
cation of the boundary of the region of periodic line deposits. 
A detailed study reveals that close to the boundary different 
patterns may often emerge for identical parameter values, de- 
pending on the initial condition. In other words, many of the 
transitions are hysteretic: If one starts inside the region (h) of 
Fig.[9]with a simulation that gives regular lines and then moves 
by small parameter increments closer to the boundary one can 
detect at which parameter values the deposit turns into a layer 
of constant thickness, and in this way define a boundary of 
region (h). However, alternatively one may start in the region 
outside (h), where after some initial transient one obtains a flat 
deposit layer, and then slowly move towards the boundary of 
(h). Detecting at which parameter value the flat deposit turns 
into a regular line pattern, one finds that this transition occurs 
at a point inside the region (h), i.e. there is a small hysteretic 
region where both patterns may occur. Using this technique 
we numerically obtain the width of the hysteresis region on 
the cuts through the (i2o,0o) parameter plane that were dis- 
cussed in section |3.3| and indicated in Fig. [9] The results for 
the width of the hysteresis region are listed in Table [TJ 

From our investigations of the transitions in the (X2o,0o) 
plane, we find four typical scenarios for the onset of the for- 
mation of periodic line patterns: 

(i) At low 0o and small or intermediate Qo the spatial pe- 
riod of the lines diverges while the line amplitude first 
slowly increases and then converges to a finite value (see 
Fig. [Tol l. This indicates the occurrence of an infinite pe- 
riod bifurcation, that could be either a SNIPER (Saddle 
Node Infinite PERiod) bifurcation or a homoclinic bi- 
furcation.^ The fact that for Q = 4.64xl0~ 7 we do 
not see any hysteresis (see Table [TJ points towards a 
SNIPER bifurcation, as a homoclinic bifurcation often 
involves some hysteresis between a stable steady or sta- 
tionary state (in the present case, the deposition of a flat 



layer) and a stable time-periodic state (here, the deposi- 
tion of line patterns). Note, however, that we are not able 
to come sufficiently close to the bifurcation point to test 
whether the typical power law relation between the pe- 
riod and the distance to the bifurcation point also holds 
here. A similar behaviour is found at high Qo and fixed 
00 = 0.31 (Fig. 12 1. We discuss further at the end of this 



section how this finding compares to results for depin- 
ning transitions in other soft matter systems. 

(ii) For high evaporation rates Qq and high 0o the line ampli- 
tude decreases with increasing Qo as the boundary of re- 
gion (h) is approached, before suddenly jumping to zero 
(see Fig. 12 1. At the same time, the line period tends 
to a finite value. Outside region (h), the pattern ceases 
to be periodic and the simulations show an initial tran- 
sient deposition of a finite number of lines, followed by 
a flat layer. Since the hysteresis here is rather large (see 
Table [TJ, we conclude that the transition most likely cor- 
responds to a subcritical Hopf bifurcation. 

(iii) For low Qq and all but the very small values of 0o, as 
the boundary is approached the line amplitude first de- 
creases (with increasing rate) with decreasing Qo, before 
suddenly jumping to zero (see Fig. [T2|>. At the same time, 
the line period tends to a finite value. Just outside this 
boundary to region (h) the simulations exhibit a transient 
deposition of a finite number of lines, followed by a flat 
layer. Since there is some hysteresis (see Table[JJ), at first 
sight the transition seems to correspond to a subcritical 
Hopf bifurcation, and therefore seems to be rather similar 
to the case described in the previous point (ii). However, 
at intermediate 0o, close to the boundary of (but within) 
region (h) one sees signs of a period doubling and there 
is also the band of intermittent patterns of type (g) just 



outside of region (h) [cf. Fig. 13 (left column, first panel) 
and Fig.[8jg) and (e), respectively]. This indicates that 
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the transition might involve several complex eigenvalues 
and be related to the intermittency and period-doubling 
route to chaos J^J Whatever is happening at the boundary, 
it is certainly more complex than case (ii). 

(iv) The last case we mention is a hypothetical transition sce- 
nario that we did not see but that we believe is likely to 
occur. Our study is based on a model that can be com- 
putationally expensive, especially for values of (X2o,0o) 
that are close to the onset of the deposition of periodic 
lines. Table [T] shows that the width of the hysteresis re- 
gion Aabs varies along the boundary. In particular, at 
large concentrations Agbs becomes very small. We be- 
lieve it is likely that there may exist a small interval along 
the boundary where there is no hysteresis, i.e., A a ^ s = 0, 
and the line amplitude gradually decreases to zero while 
the line period approaches some non-zero value, and the 
deposition pattern becomes a small amplitude harmonic 
modulation. Such a transition would correspond to a su- 
percritical Hopf bifurcation. 

Note, that the present study only considers one-dimensional 
deposition patterns. Two dimensional deposition patterns are 
beyond our scope. We expect the full two-dimensional be- 
haviour to be very rich. In particular, at small evaporation 
rates one should expect the evaporative dewetting front to be 
transversally unstable, even in the situation without solute.^ 
We expect this to also occur in the case with solute, close to 
the transition discussed in points (i) and (iii). This argument is 
bolstered by the observation that the particular transition de- 
scribed in (i) involves deposits that are nearly homoclinic in 
space. There exist generic results^ that show that patterns 
near homoclinic solutions are prone to instabilities. 

Before we move on to discuss in the next section the in- 
fluence of system control parameters besides the solute con- 
centration <pQ and the evaporation number Qq, we make a few 
comments to put our findings into a wider context. It is im- 
portant to understand that the observed transitions from sta- 
tionary front motion to the deposition of lines may be seen as 
depinning transitions in the frame moving with the mean front 
speed: When a front of constant speed deposits a flat layer, in 
the comoving frame the concentration profile is steady. Then, 
one may say that the concentration profile is pinned to the 
moving front as it does not move relative to it. However, at 
the transition to depositing a periodic pattern, the concentra- 
tion profile starts to stay behind the moving front, and one may 
say it depins from the front. Note, however, that after depin- 
ning, the concentration profile does not move relative to the 
front as a whole. Instead, only a part of it (the jammed part) 
starts to move relative to the front, resulting in the deposition 
of a line. This process then repeats periodically. 

From this observation it becomes clear why the transition 
scenarios (i) to (iv) described above are analogous to similar 



scenarios found in studies of depinning in other driven soft 
matter systems. To illustrate how universal such transitions 
are, we mention three systems: First, drops of simple non- 
volatile liquids that sit on heterogeneous substrates and de- 
pin from the heterogeneities under the influence of external 
driving forces. Depending on the particular setting and pa- 
rameter regime, one may observe SNIPER, homoclinic and 
super- or subcritical Hopf bifurcations.^^MSI There, however, 
the entire depinning drop or ridge slides along the heteroge- 
neous substrate in a periodic manner. A second system con- 
sists of clusters of interacting colloidal particles that shuttle 
under the influe nce of external forces through a heterogeneous 
nanoporeJ 89 ! 90 ! Under weak dc driving, the peak in the parti- 
cle density distribution is pinned by the heterogeneities of the 
pores. However, depending on driving force and the attraction 
between the colloids, depinning transitions via Hopf and ho- 
moclinic bifurcations occur resulting in time periodic fluxes. 
The third system is much closer related to the one studied here, 
as it also concerns patterns that are produced at a three-phase 
contact line: If a Langmuir-Blodgett surfactant monolayer is 
transferred from a bath onto a solid substrate, one may ob- 
serve stripe patterns in the deposit resulting from substrate- 
mediated condensation.!^"^ Time simulations of a dynami- 
cal model for this system show the occurrence of patterns of 
stripes parallel to the contact line (and other patterns too)J^ 
In a reduced model the depinning transitions from a steady 
(pinned) concentration profile to the time periodic (depinned) 
state may occur via a homoclinic bifurcation or a subcritical 
Hopf bifurcation.!^ 

The comparison with these different soft matter systems 
that show depinning, corroborates the picture we have given 
above in points (i) to (iv). However, a systematic analysis of 
the various transitions related to the onset of the deposition of 
line patterns in the present system requires a simplified model 
to be developed. 

In the following sections we briefly describe how varying 
other physical aspects of the system influence the behaviour 
of the system. In Sec. |3.5| we examine the effect of varying the 
diffusivity of the solute in the solvent. In Sec. |3.6| we discuss 
the influence of varying the wettability and also of varying the 
chemical potential of the solvent in the ambient vapour. In 



Sec. 3.7 we discuss effects of solvent rheology. 



3.5 Influence of diffusion 

The influence of diffusion is quantified in our model by the in- 
verse Peclet number Pe _1 . Up to this point we have discussed 
the basic mechanism as being based on a subtle balance of 
convective and evaporative motion. In the results presented so 
far, the solute diffusivity was low and contributed little to the 
overall transport. 



Fig. 14 presents results for several measures which charac- 
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Fig. 14 (Color online) Measures characterising the regular line 

patterns as a function of Pe~ 1 for fixed parameters: 82 * (purple 

line), the standard parameter values; o (blue line), 

Qq = 1.667 x 10 _7 ,^o = 0.41; A (green line), 

Qq = 4.64 x 1(T 7 , O = 0.3016; (red line), 

Q Q = 5.663 x 10~ 6 , O = 0.41; and □ (light blue line), 

£2q — 4.64 x 10~ 7 , 0o = 0.498. The purple line corresponds to states 

in the central part of the region (h) on Fig. [5] whereas the other lines 

correspond to points close to its boundary. Selected corresponding 

line patterns are shown in Figs . 1 1 5 1 and| 1 6| Note, that for A (green) 

the period P monotonically increases to very large values with 

decreasing Pe _1 (up to P « 30000 at Pe" 1 = 0). 



terise the periodic line patterns, as a function of the inverse 
Peclet number Pe _1 , which is a ratio of the time scales for 
convection and diffusion. Results are shown for five sets of 
values of the parameters (Qq, 0o), corresponding to locations 
within the region of periodic lines [region (h) in Fig. [9): sym- 
bols * indicate our standard values of (X2o,0o) an d the other 
sets correspond to locations close to the boundary of region 
(h). Examples of the corresponding line patterns are displayed 



in Figs. 15 and 16 All the other parameters, %, M, and v are 
equal to the standard values (see the end of Section[2]). 

Around and below the standard value of Pe _1 = 0.0003, 
the deposition pattern is almost independent of the value of 
Pe _1 . Decreasing Pe _1 to zero has almost no effect on the 
size of region (h), and the only effect is that the deposit pat- 
terns becomes slightly sharper. However, on increasing Pe~ 1 
above 0.0003, the size of the region (h) (Fig. [9]) starts to shrink 
considerably, until it vanishes entirely as the effects of solute 
diffusion increase (roughly, when Pe~' > 0.003). The shape 
of the lines changes monotonically and they become more si- 
nusoidal in shape: the amplitude, period and skewness all be- 
come smaller with increasing Pe _1 , while the relative width, 
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Fig. 15 (Color online) Morphology changes of the line pattern as 
Pe~' is varied, corresponding to Fig.[l4] The left column shows 
profiles corresponding to the blue line with o symbols and the right 
column to the purple line with * symbols in Fig. [14] From top to 
bottom, the panels in both columns correspond to 
Pe" 1 = 3x 10" 5 ,10" 4 ,1.7x 10~ 4 ,3x 10~ 4 .5.6x 10" 4 and 0.001. 



K = 2(7 1 P increases (see Fig. 14 1. We remind the reader that 



highly peaked or separated lines have small K values, whereas 
a deposit profile that resembles a harmonic wave has a large 
K value. When the diffusive mobility of the solute is large 
the effect of diffusion counteracts the solute build up due to 
evaporation. The transition between the deposition of peri- 
odic lines and a flat layer occurs in between Pe~' = 0.001 and 
Pe~' = 0.003 for the standard values of (i2o, (j>o) and at some- 
what smaller Pe~' for parameter values close to the boundary 
of region (h), as the region itself is shrinking. Such a large dif- 
fusivity is unlikely for large nanoparticles in suspension but 
might occur for very small particles or for small molecules 
in solution. Since the line amplitude approaches zero at this 
transition point, whilst the period remains finite and the line 
profiles become nearly harmonic, we suspect that this transi- 



tion corresponds to a Hopf bifurcation (cf. section 3.4 1 but we 
did not study this transition in as much detail as the transitions 
discussed in Sec. 13.41 



In the lower right panel of Fig. 14 we display results ob- 
tained for the skewness, 73. As Pe _1 is increased, the curves 
marked by symbols o and □ show at first an increase of "fy be- 
fore it decreases again just before the pattern vanishes. The 
reason for this non-monotonic behaviour can be seen in the 



left column of Fig. 15 that displays deposit profiles which cor- 
respond to the o curve in Fig. 14 There, going from the first 



to the fourth panel from the top we see new secondary lines 
emerge out of the tail of the primary lines, i.e., decreasing 
the effects of diffusion, leads to a period doubling. The period 
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Fig. 16 (Color online) Morphology changes of the line pattern as 
1 is varied, corresponding to Fig. 
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The left column shows 



profiles corresponding to the green line with A symbols and the 
right column to the red line with symbols in Fig.fW] From top to 
bottom, the panels in both columns correspond to 
Pe" 1 = 3 x 1(T 5 , 1CT 4 ,3 x 1(T 4 and 0.001. 



doubling is not visible in the dependence of the period on Pe 



(see Fig. 14 1 because around these parameter values the pat- 
tern is not very regular and the secondary lines appear through 
a smooth transition. These irregular patterns seem to be analo- 
gous to the intermittent patterns [Fig.|HJg)] or the depression- 
line patterns [Fig. [8]j)] but here these intermediate patterns 
appear in a somewhat wider band of parameter values close to 
the border of region (h) in the phase plane shown on Fig. [9] 
The period of these irregular patterns P s» ishift/2, which in- 
dicates that the shift length (see Appendix [A} interacts with 
the period of the pattern. In general, all such irregular patterns 
appear only close to the boundary of region (h) (Fig. [9} of pe- 
riodic lines, where Qq is small and (j>o is moderate to high. In 
other words, where patterns have large positive 73 (heavy right 
tail), we observe the 'chaotic deposition' described in Sec. 
and we expect transversal instability effects. In Appendix 



3.4 



A.3 



we discuss how robust the results we obtain are, in relation to 
the numerical methods we use to solve our model equations. 



The profiles in Figs. 15 and 16 also show that as the value of 
1 is increased, the parameter regions where lines of large 



Pe 



positive skewness (long tail towards front) are found (smaller 
Qo and moderate and high 0o) become replaced by lines that 
are nearly symmetric (small skewness) or point the tail away 
from the front (negative skewness). In addition, the period of 
the lines decreases. 

Finally, we mention the importance of diffusion in ex- 
periments and applications of such evaporation driven self- 
organisation processes. The relative importance of diffusion 



depends on both the mobility of the solute and the thic kness 
of the solution layer. For instance, for very thin filmsEIISl 
one sees from Eqs. Q and d3), that the diffusive mobility 
(2d 06 h dominates the convective mobility Q c oc /j 3 ; see also 
the discussion by Vancea et al. where the dynamics in ultra- 
thin postcursor films left behind mesoscopic dewetting fronts 
of nanoparticle suspensions^ is modelled by either kinetic 
M onte C arlo model s 92 " 94 ' or dynamical density functional the- 
oryl 95 ! 96 ! that only incorporate diffusive transport because in 
this situation diffusion is the dominant process. The present 
thin film model captures the competition of convective and dif- 
fusive transport relevant for larger film thicknesses where the 
key influence is the mobility of the solute itself. Xu et al.&i 
found that as the size of nanoparticles is decreased, one ob- 
tains spoke-like structures and irregular patterns instead of the 
regular rings seen for larger nanoparticles. This indicates that 
an increase of the influence of diffusion brings the system 
closer to the border of the region of periodic line deposition, 
where patterns become irregular and unstable in accordance 
with our present findings. 

In summary, a low value for the solute diffusivity means 
diffusion does not influence the deposition patterns and diffu- 
sion may actually be neglected in the model. However, fast 
diffusion is able to counteract evaporation and effectively sup- 
presses the occurrence of deposition patterns. 

3.6 Influence of wettability and chemical potential 

After having discussed the influence of solute diffusion in the 
previous section, here we briefly consider the influence of wet- 
tability and the chemical potential of the solvent vapour. We 
start with the influence of the wettability that is quantified in 
our model by the parameter %, contained in the polar con- 
tribution to the disjoining pressure Tl = l/h^ — exp(— %K). 
The influence of changes in % on n can be appreciated by 
inspecting Fig. [3] A decrease [increase] from our standard 
value % = \ .085 results (at fixed chemical potential M) in a 
decrease [increase] of the precursor film thickness and of the 
thickness of the bulk film; an increase [decrease] of the thick- 
ness where the front profile has its inflection point; a decrease 
[increase] of the wetting energy at the precursor film height, 
and in consequence an increase [decrease] of the slope at the 
inflection point, that is our measure of a 'nonequilibrium con- 
tact angle' . It also leads to larger [smaller] energy difference 
between the two stable heights, h\ and h%, so the evaporation 
will be stronger [weaker] . 

Corresponding results for the deposition patterns are dis- 



played in Fig. 17 (line characteristics) and Fig. 18 (selected 



deposition profiles). The line amplitude and period decrease 
monotonically with increasing %, while the skewness first in- 
creases and then drops in value. The period, amplitude and 
relative width behave in a similar manner as when increasing 
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Fig. 17 (Color online) Measures characterising the regular lines as a 
function of ^.^HThe blue lines with o symbols have fixed 
Qq = 7 X 10~ 7 , and the standard concentration 0o = 0.41, whilst the 
green lines with □ symbols have the (standard) evaporation rate 
£2g = 4.64 x 10~ 7 and concentration 0y = 0.31. All other parameter 
values are the standard ones, defined at the end of Sec. [2] The blue 
lines corresponds to configurations in the central part of the region 
(h) on Fig. [9] Selected line patterns corresponding to these are 
shown in Fig.[T8] Note that the period P seems to decrease 
exponentially with increasing % but we also note that for % smaller 
than the standard value, we observe stable double lines. 



00, see Sec. 3.3 This indicates that the transition towards a flat 



scenario (ii) or (iv) discussed in Section 3.4 On decreasing %, 
we find that the lines become more anharmonic and sharper 
and after a period doubling become deposited in pairs. A fur- 
ther decrease results in groups of multiple lines and in general 
to the deposition of a more irregular pattern, i.e. the transition 
towards a flat deposit at small % is most likely via scenario 
(iii) introduced in Section 3.4 Some corresponding profiles 



are displayed in the left hand column of Fig. [T8] 

The main features of an evaporative dewetting front are in- 
fluenced not only by x, but also by the dimensionless chemical 
potential M (cf. Fig. [3}. A decrease in M at fixed % results in 
a decrease of the precursor film thickness and also of the up- 
per film thickness; an increase of the thickness where the front 
profile has its inflection point; a small increase of the wetting 
energy at precursor film height but also in an increase of the 
energy difference between the two stable heights, hi and h%; 
and, as a consequence an increase of the contact angle. A 
further decrease in M leads to faster evaporation in the front 
region, i.e., a larger front velocity, and so the capillary ridge 
decreases and the 'nonequilibrium contact angle' decreases. 
Many of these features can already be seen in the case without 
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Fig. 18 (Color online) Morphology changes as % is varied, 
corresponding to the results in Fig.[T7] Left column: the profiles 
correspond to the blue lines with o symbols in Fig.|17| from top to 
bottom the results are for % = 1.06, 1.062, 1.065, 1.075 and 1.095. 
Right column: corresponding to the green lines with □ symbols in 
Fig.|17| from top to bottom the results are for 
X = 1.062, 1.065, 1.075, 1.085 and 1.095. The other parameters 
have the standard values. 



deposit at large % is most likely via a Hopf bifurcation, i.e. via solute 121 



Results obtained for varying M are presented in Fig. 19 A 
decrease in M, i.e. increase in |M|, results in an increase of the 
line amplitude and period. We do not consider M < —0.0033 
as then no second stable film height exists, a case not covered 
by our numerical set-up. Increasing M, the line amplitude de- 
creases towards zero, the period approaches a constant value 
and the relative width and skewness of the lines both increase. 
This is similar to the case of increasing % (compare the two 
columns of Fig. 17 1. We see this as an indication of a Hopf 
bifurcation as the most likely transition mechanism, i.e. sce- 
nario (ii) or (iv) of Section [3~4] If this is correct, the skewness 
should decrease strongly over a small range of values of M. 
Note finally, that it is difficult to separate the influences of 
the parameters % and M as they both influence the stable film 
heights, the contact angle, and whether a capillary ridge exists 
or not. 



3.7 Influence of solution rheology 

In this final results section we discuss the influence of solution 
rheology - i.e. we present results for varying the exponent v 
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Fig. 19 (Color online) Left column: Measures characterising the 
regular line pattern as M is varied. 82 The upper panel shows both 
the amplitude (blue line with o symbols) and the period (green line 
with □). The Lower panel shows the relative width (blue line with o) 
and skewness (green line with □). All the parameters except M have 
the standard values as defined at the end of Sec. [2] Right column: 
Corresponding profiles for ifcf = -(3.3,3.2,3.1,2.9,2.8) x 1(T 3 , 
from top to bottom. 



in the Krieger-Dougherty law (fT6]). In the thin film literature 
the value v = 2 is often usedJ^MTand our standard configura- 
tion employs the value V = 1 .575 corresponding to spherical 
colloidal particles as discussed above following Eq. Q. How- 
ever, some authors suggest that for some systems, much lower 
values are required, in particular, in the context of jamming 
transitions for attractive colloidal particles. There, values as 
low as v = 0.13 are mentioned.^ 

Here we consider moderate variations in the interval 1 < 
V < 2, and investigate how the deposition patterns change for 
five selected points in the parameter plane (£2o,0o)- Three of 
the configurations correspond to the crossing points of the four 
straight dashed lines in Fig.|9]that lie within the region (h) of 
periodic lines [one is the standard case described at the end of 
Sec. [2] and the other two lie close to the boundary of region 
(h)]. The other two selected points also lie close to the bound- 
ary of the region (h). Results for measures characterising the 
deposition patterns are displayed in Fig. [20] and selected cor- 



responding deposition profiles are displayed in Figs. 21 and 

E2] 

First, we consider the case of the standard set of parameter 
values, that is in the central part of the region (h) in Fig. [9] 
The measures are shown with the symbols o (blue curve) in 
Fig. [20] and selected corresponding profiles are displayed in 
the left column of Fig. 




Fig. 20 (Color online) Measures characterising the line pattern as a 
function of the viscosity exponent v.^The blue curve with o 
symbols is for fixed standard £2q = 4.64 x 10~ 7 and standard 
tj){) = 0.41; The green curve with □ is for £2q = 4.64 x 10~ 6 and 
standard 0o; The red curve with is for the standard value of Qq 
and for 0y = 0.31; The light blue line with A is for the standard 
value of Qq and for tj>o = 0.498; The purple line with * is for 
£2q = 1.667 x 10 -7 and standard <j>Q. All the other parameter values 
are standard, as defined at the end of Sec. [2] The blue line 
corresponds to configurations in the central part of the region (h) on 
Fig. [5] whereas the other lines correspond to points close to its 
boundary. The corresponding line patterns are shown in Figs. |21| and 
1331 



ative width decreases and the skewness departs from zero in- 
dicating that the lines become more peaked and anharmonic at 
lower v. Note that simulations tend to be more computation- 
ally demanding at lower v and e.g., the point at v = 1.0 has 
a larger error, because effects of the initial transient are still 
present in the analysed data. 

The other selected points, located close to the boundary of 
region (h) in Fig.|9]behave somewhat similarly, but with some 
differences. The symbols (red line) in Fig. [20] correspond 
to lines with a large amplitude and a long period; a selection 
of profiles for this case are displayed in the right column of 
Fig. 21 Interestingly, the period changes nonmonotonically, 
while the skewness decreases as v is decreased. Profiles that 



correspond to the symbols □ (green line) in Fig. [20] are dis- 
played in the left column of Fig. 22 The lines have negative 



21 As v decreases from the value 



V = 2, the amplitude of the deposition lines increases, the rel- 



skewness and vary with v in a manner similar to the results 
for the standard set of parameter values. The symbols A (light 
blue line) and * (purple line) correspond to the points close to 
the boundary of region (h) where period doubling transitions 
are observed and results become quite sensitive to the details 
of the numerical solution method. The computations are very 
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Fig. 21 (Color online) Morphology changes for various values of 
the exponent v, corresponding to Fig. [20] The left column 
corresponds to the blue lines with o symbols and the right column to 
the red lines with in Fig.[20] From top to bottom panels in each 
column correspond to v = 1, 1.25, 1.5, 1.575, 1.75 and 2. 



sensitive for decreased v values for the case when 0o = 0.498, 
with standard Qq, which correspond to the symbols A (light 



blue line) in Fig. 20 with corresponding profiles displayed in 



the right column of Fig. 22 There we see patterns that switch 



between a depression-line pair and a flat layer, and irregular 
intermittent patterns. 

In general, the deposition of regular lines is robust with re- 
spect to a change in v and we observe them over a broad in- 
terval of parameter values. This is to be expected from the 
experimental results, where periodic lines were seen for var- 
ious (non)spherical nanoparticles,^ colloids^ and polymeric 
solutes. 32 The general trend is for the regular line patterns to 
become more sharply pronounced as the value of v becomes 
smaller, i.e., the line amplitude increases and their relative 
width decreases. 

We do not study in detail the onset of periodic deposits 
when varying v. For very small values of v the simulations 
are computationally too demanding, since the lines are more 
sharply peaked. For some configurations close to the bound- 
ary of region (h) in Fig. [9] we observe irregular deposits and a 
transition to a flat deposit indicating that region (h) may shrink 
as v is decreased - cf. Fig. 22 As the value of v is increased, 



the line amplitude decreases, and we expect that for values of 
V larger than 2, the region of periodic deposits will also van- 
ish. 
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Fig. 22 (Color online) Morphology changes for various values of 
the exponent v, corresponding to Fig. [20] The left column 
corresponds to the green line with □ symbols and the right column 
to the light blue line with A in Fig. [20] From top to bottom panels in 
each column correspond to v = 1, 1.25, 1.5, 1.575 and 1.75. 



4 Discussion and conclusion 

We have presented a thin film model for the close-to- 
equilibrium self-organised deposition of material onto a 
smooth flat solid surface from a receding three-phase con- 
tact line of a polymer solution or nanoparticle suspension. 
The model consists of a pair of highly non-linear coupled 
long-wave evolution equations for the film thickness h and 
the effective solute layer height h v = h(j), where is the 
height-averaged scaled solute concentration. The evolution 
of the film is driven by a number of terms in the equations 
which describe several physical effects, including: (i) capil- 
larity through a Laplace (or curvature) pressure, (ii) wettabil- 
ity through a Derjaguin (or disjoining) pressure, (iii) evapora- 
tion due to a difference in the local solvent chemical poten- 
tial and that of the vapour, and (iv) forces due to gradients in 
the solute concentration. The transport processes that are in- 
volved in the dynamics correspond to convective and diffusive 
transport (both give a conserved dynamics), and evaporation 
(a non-conserved dynamics). An important ingredient in the 
model is the rheological property of the solution/suspension, 
that leads to an arrest of the convective motion at some crit- 
ical solute concentration, e.g., at random close packing, for 
a suspension of spheres that have no net attractive forces be- 
tween them and only interact via excluded volume repulsive 
interactions. To model these effects, we have employed the 
Krieger-Dougherty power law (|6]l for the viscosity, but we ex- 
pect similar behaviour to occur for other such laws. The model 
we u se is related to some other models used in the literature 
(e.g. ,^^3, as explained in Section]^. 

Numerically solving our model equations, we have investi- 
gated the deposition of regular and irregular line patterns that 
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has been observed in numerous experiments utilising a wide 
range of materials and experimental set-ups. Assuming the 
front profiles only vary in one spatial direction, we have found 
that regular line patterns are formed over an extended region 
of the parameter space. Other one-dimensional patterns that 
we have encountered, include the transient deposition of a sin- 
gle or a finite number of lines, periodic arrays of double lines, 
a periodic switching between depression-line pairs and a fiat 
layer, and irregular intermittent line patterns. 

We believe that the model explains one of the basic mech- 
anisms for the formation of regular line patterns. They re- 
sult from a self-organised cycle of deposition-caused pinning- 
depinning events t hat is e xperimentally often described as a 
'stick-slip' motionJ 22 ! 26 ! 30 The stick-slip motion is caused by 
the highly nonlinear rheology: Evaporation leads to a rapid 
increase of the solute concentration in the contact line region 
resulting in a strong increase of the viscosity. This in turn 
eventually leads to arrest of the convective motion of the re- 
ceding front. However, the front is not entirely stuck, it still 
moves due to evaporation, albeit at a much slower speed. Dur- 
ing this phase more material is deposited resulting in a line 
deposit. After a sufficient amount of material is deposited, the 
concentration in the contact line region decreases, the front 
depins from the line deposit, and starts moving with much 
greater speed. This appears to be a stick-slip motion. Thus, 
the self-organisation of the patterned deposit results from a 
subtle interplay of all three of the transport processes (con- 
vection, diffusion and evaporation) and the stick-slip motion 
results from the large difference in the timescales of convec- 
tive and evaporative de wetting. 

After establishing this basic mechanism, we have per- 
formed a detailed study to investigate the influence of solute 
concentration, evaporation rate, diffusion, wettability, chem- 
ical potential, and solution rheology on the patterns that are 
formed. In general, line patterns emerge at intermediate values 
of both the evaporation rate and the solute concentration, indi- 
cating that both are important quantities in the deposition. Fur- 
thermore, line patterns are suppressed when the solute diffu- 
sivity is sufficiently high, because when this occurs, the mov- 
ing front is not able to collect the solute. In the opposite limit, 
when solute diffusion is set to zero, we find that this does not 
significantly affect the line patterns. In this paper, we have 
probed the influence of wettability solely through the param- 
eter x m the short-range part of the Derjaguin pressure. Our 
results indicate that an increase in the wettability may lead to 
a suppression of the patterns. Note, however, that our inves- 
tigations of this aspect are not exhaustive. We have not fully 
studied this aspect here, because our present numerical set-up 
does not allow for the study of a completely wetting solution. 
Such a study will be pursued in the future, employing a mod- 
ified set-up for the case of an active geometry. The influence 
of the solute on the system rheology is quite notable but not 



easy to categorise: decreasing the exponent v in the power law 
relation between the viscosity (and diffusivity) and the solute 
concentration, increases the amplitude of the lines for all the 
values of the other parameters that we have tested. However, 
the influence on the period of the line pattern depends on the 
other parameters: it can increase, decrease or even behave in 
a non-monotonic manner with decreasing v. 

In gene ral, thes e findings are in agreement with experimen- 
tal results.SSSHHH] However, to our knowledge, no experi- 
mental study has systematically mapped out the regions where 
the various line patterns exist (analogue to our Fig. [9]), or pre- 
sented a detailed (quantitative) analysis of how the line pat- 
tern characteristics depend on the system control parameters. 
Therefore, before we make a comparison of the typical de- 
pendencies found in the experiments and our model, we first 
discuss the typical length and time scales in the deposition 
process. 

The scales we have used [cf. Eq. ( [TT| ] for the time, x- 
coordinate, and film thickness are X = 3rjoY/8\S p \ 2 , a = 
{8y/\S p \) l l 2 , and 8 = {A/6n\S p \) l l i , respectively. For our 
standard parameter values (see Sec. |2), we find T = 2.3 x 
10~ 8 s, a = 2.9nm, and 8 = 0.5nm, i.e., the line pattern 
in Fig. [5] has a period of A sa 3/im and the mean deposition 
velocity is v 1/im/s. For Fig. [6] we have X ss 3/zm and 
v « 0. 1 jUm/s, whereas the periods of the patterns displayed in 



Fig. 10 range from about 1 .7/im to about 52/im. These values 
are similar to those of the experiments by Xu et aZ.^S. Their 
Fig. 1 shows line patterns with 3-5/zm periods. Fig. 2 of Xu 
et al. Ogives line periods between 4 and 20/im, a range sim- 
ilar to results by Hong et alf^l. Some characteristics that we 
find are rather similar to experimental findings by Bodiguel 
et al.^: The contact line velocity is similar to the 0.8/xm/s 
imposed in their Fig. 1 ; the spatial periodicity of the pattern in 
their Fig. 2 is 210/xm and the contact line speed is 6.2/im/s; 
In general, their typical imposed velocities are 1-10/im/s. 

It is hard to compare with experimental results in the liter- 
ature how pattern properties typically depend on the control 
parameters, as most works do nearly not discuss these. Some 
authors discuss dependencies that we are not able to measure 
in our simulations such as, e.g., the pinning forced A de- 
crease in the ring spacing with decreasing distance from the 
probe center is sometimes described in studies of evapora- 
tive dewetting with the meniscus technique (see Figs. 1 and 
3 of Kwon et al.^, Figs. 2 and 3 of Hong et al.^1, Lin and 
Granick 97 ), and Fig. 2 of Xu et al™. In all mentioned cases 
this may be interpreted as observing a decrease in the period 
with increasing solute concentration and/or decreasing evap- 
oration rate. Both of these possibilities are in accord with 
our results in Figs. 10 and 12 respectively. However, a de- 
crease in the period with decreasing concentration is found by 
Xu et al. 25 , a result also reported by Xu et al. 22 when com- 
paring evaporating polymer solutions of different initial con- 
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centrations. Adachi et alM reports a non-monotonic depen- 
dence of the line density on the initial solute concentration. 
Decreasing the concentration leads to the line density slowly 
increasing before abruptly decreasing at the edge of the pat- 
terning region. This is not unlike our findings - see Fig.[lO] A 
quantitative comparison is difficult, as they report observing 
about 1-8 rings produced with a small evaporating droplet. 
Related problems obstruct a detailed comparison with many 
of the experiments performed with circular drops, which have 
radial geometries, that always imply a drift in several system 
parameters in the radia l direction. The planar geometry used 
by Bodiguel et a/J 26 j 98 l allows for a more natural comparison. 
However, there detailed results are given for the dependence 
of the pinning force on experimental parameters whereas no 
such results are given for the period of the line patterns. What 
is remarkable, is the agreement for the morphology of indi- 
vidual lines: Bodiguel et al.^ shows a rather asymmetric line 
profile with a negative skewness, i.e., the long tail points away 
from the receding contact line, similar to the profiles displayed 
in last two panels of the left column of our Fig. 13 Such 



agreement also exists with results for some unpublished ex- 
periments with nanoparticle suspensions 83 where line profiles 
having a negative skewness and also a positive skewness were 
found, depending on the solvent that was used, i.e., the long 
tail can point either away or towards the receding contact line, 
respectively, similar to the transition we see that is illustrated 
in Figs.[l2]and[T3] 

Within the limits of our numerical approach, particular care 
was taken in the analysis of the onset of the formation of lines 
when changing the various system control parameters. We 
have observed three main types of transition at the border of 
the region of periodic lines: (i) the line period diverges while 
the line amplitude converges to a finite value, (ii) the line pe- 
riod doubles, followed by a region of intermittent patterns, (iii) 
the period approaches a finite value while the amplitude first 
decreases before suddenly jumping to zero. We have also pro- 
posed a hypothetical fourth scenario - (iv) a case where the 
line period approaches a finite value while the line amplitude 
approaches zero. To relate these findings to dynamical sys- 
tems or bifurcation theory, it is helpful to consider a stationary 
moving front that deposits a homogeneous layer as a steady 
state in a comoving reference frame. At the onset of the line 
formation this steady state becomes unstable and gives way 
to a time-periodic solution (in the frame that moves with the 
mean speed of the front). In this context, the above men- 
tioned transitions can be identified as (i) an infinite period bi- 
furcation (either homoclinic or Saddle-Node Infinite PERiod 
(SNIPER)), (iii) a subcritical Hopf, and (iv) a supercritical 
Hopf bifurcation, respectively. Case (ii) is more complicated 
and resembles the period doubling and intermittency route to 
chaos when approaching the transition from the side of the 
periodic deposits. Approaching the transition from the side of 



the flat deposit, one finds that the length of the transient seems 
to diverge, which is similar to the transition from an excitable 
medium to one that exhibits sustained oscillations. In particu- 
lar, the homoclinic bifurcation and the subcritical Hopf bifur- 
cation have recently been identified to be responsible for the 
transition to stripe-like deposition patterns in the Langmuir- 
Blodgett transfer of a surfactant layer onto a moving plate, 
where patterning is due to substrate-mediated condensation of 
the surfactant.^EH 

To put these transitions in a wider context, we would like 
to point out that the onset of stripe formation can also be seen 
as a depinning transition: When a homogeneous layer is de- 
posited from a stationary moving front, the concentration pro- 
file of the solute is pinned to this front. However, this is not the 
case when lines are deposited. In this case, the lines are de- 
pinned from the average-velocity comoving frame - i.e. they 
stay behind the moving front. Note however that many details 
of the present transition are unlike the depinning via a homo- 
clinic or SNIPER bifurcation described by Beltrame et al.^ 
who study how a drop depins from a substrate heterogeneity 
under the influence of an external driving force, leading to the 
entire drop sliding along the substrate. Here, after depinning, 
the concentration profile at the moving front does not move in 
its entirety away from the front. Instead, only a part of it (a 
line) starts to move - a process that is then repeated periodi- 
cally. Note that for other parameter values for the driven drop 
system, the depinning transition may also be due to a Hopf 
or a SNIPER bifurcation.^ The present system appears to be 
more complicated, as the three main types of transition can 
be altered by interactions with period-doubling bifurcations 
and/or intermittent behaviour. We believe it is necessary for 
reduced models to be developed in order to investigate these 
issues in detail. Although our understanding of the onset of 
the patterning is still not complete, we hope that this aspect 
of our investigation may prove useful in the classification of 
future experimental results for the onset of the formation of 
deposition patterns. 

In the present work we have limited our attention to a one- 
dimensional geometry, i.e., we have studied deposition from a 
receding straight front with an imposed transversal stability of 
the front. A fully two-dimensional treatment is worth pursu- 
ing (but computationally not feasible using the numerical ap- 
proach taken here - see appendix[A| and highly relevant, since 
it is well known that receding contact lines may be transver- 
sally unstable, in particular, if the recedi ng liquid is a suspen- 
sion or solution with a volatile solvent.D- 2 ! 1 3 I 25 I 28 1 1 Q 2l Experi- 
mental results on line deposition indicate that such transversal 
instabilities are likely to occur on the left and/or upper border 
of the region of line patterns in the phase diagram we have 
presented in Fig. [9] In the experiments by Hong et a/.^l an( j 
Xu et al. ^transversal instabilities occur at very small veloc- 
ities of the receding contact line, i.e., when the evaporation 
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is slow and/or the solute concentration is high. In our model, 
we have observed that under these conditions the deposition 
process shows period doubling and results become very sen- 
sitive to numerical details, which suggests the system may be 
transversally unstable, in agreement with these experiments. 
In the system studied by^M different types of transversal insta- 
bilities occur at low and high solute concentrations. However, 
this system represents an example of an active geometry where 
the sliding velocity of an upper plate controls the patterns that 

1981 

are formed. Karthaus et cuM® states that a Plateau-Rayleigh 
instability of the rim at the moving front is responsible for 
the fingering instability, an observation which is in line with 
arguments made for receding fronts of non-volatile^^ and 
volatile 75 pure liquids. However, the name "Plateau-Rayleigh 
instability" might be misleading in the case of a moving front 
of a solution or suspension since (i) the motion of the liq- 
uid rim itself may change the character of the instability,!^ 
(ii) p artia l wettability may trigger an instability of a receding 
front,E2H an( j 5 most importantly, (iii) since the solute concen- 
tration may change along the front, new classes of instabilities 
are possible, which are related to mobility contrasts and lo- 
cal phase separation 96 or solutal Marangoni effects. These 
possibilities warrant future studies employing either geome- 
tries that allow for a more efficient numerical approach and/or 
the development of reduced models in a similar spirit to the 
model for Langmuir-Blodgett transfer of a surfactant layer 
developed by Kopf et al)^. A thin film model for this pro- 
cess predicts a transition from stripes that are perpendicular 
to the receding front, to parallel, as the transfer velocity is in- 
creased.- - We expect a similar behaviour in our situation if 
our resulting front speed is taken as a control parameter. 

Our work may also be related to different approaches that 
describe other aspects of the evaporative dewetting of solu- 
tions and suspensions. In particular, the role of diffusion 
of the solute in suppressing the line deposition that we have 
described above, is very similar to its role in the suppres- 
sion of fingering instabilities that occur at receding evapo- 
rative dewetting fronts in ultrathin films of nanoparticle sus- 
pensions, that were experimentally observed in by Pauliac- 
Vaujour et al.^, and subsequently modelled by a kinetic 
Monte Carlo modefH and dynamical density functional the- 
ory JUM] These non-hydrodynamic models are not able to de- 
scribe mesoscopic hydrodynamics, i.e., the transport of mo- 
mentum, but do allow one to incorporate the various interac- 
tions of solute and solvent, i.e., of the complete thermodynam- 
ics of the system. These interactions can, if strong enough, 
result in a local phase separation of solute and solvent close to 
the moving contact line and in this manner trigger the front in- 
stability. 96 The present hydrodynamic model is not able to de- 
scribe this, but it does have the advantage that it incorporates 
basic wettability effects, i.e., the substrate-film interactions. A 
future avenue for improvements of our mesoscopic hydrody- 



namic model is to incorporate solute-solute and solute-solvent 
interaction. A related thin film model for a layer of a de- 
composing non-volatile binary mixture has recently been de- 
rived by Naraigh an d Thiff eault^^ as a long-wave approx- 
imation to model-Hj 105 " 107 ! A general approach of deriving 
such thin film evolution equations in the context of nonequi- 
librium thermodynamics (taking the form of a gradient dy- 
namic s based on a free energy functional), was recently pro- 
posed.^Mi22l This then naturally allows one to incorporate ef- 
fects like solute-dependent wettability and capillarity (includ- 
ing solutal Marangoni effects), and the dependence of evapo- 
ration on the osmotic pressure. In this way, the present model 
can readily be extended and may be employed to assess the 
influence of additional effects on the basic mechanism that we 
have describe here. This should be done for passive and active 
geometries alike. 

We acknowledge support by the EU via the ITN MULTI- 
FLOW (PITN-GA-2008-214919). 

A Numerical approach 

A.l Discretization and numerical scheme 

In order to calculate the deposit profiles h p (x,t) corresponding 
to evaporative dewetting of a thin liquid film of height h(x,t) 
of a solution or a suspension we perform time simulations of 



the model equations ( 13 1 and ( 14 1, in one spatial dimension. 
The initial film profile in the form 



h(x,0) =h\ + Q12 — h\)H(x — xq) 



(21) 



which corresponds to the sharp front, located at x = xq, that 
subsequently moves in positive-x direction, as the volatile 
component evaporates. The heights hi and hi denote the 
precursor and upper stable film thicknesses, respectively, see 
Fig. [5J and H(x) is the Heaviside step function. The initial 
concentration profile is constant, <j>(x,0) = h(x,0)/h p (x,0) = 
00- Furthermore, to approximate the semi-infinite drying thin 
film on a finite computational domain, particular boundary 
conditions (BCs) are employed, as discussed later in the Sec- 
tion A. 2 The computational domain x <G [0,1] is discretized 



employing an inhomogeneous grid adapted to resolve the de- 
tails in our unknown fields h(x.t) and h p (x,t). For time step- 
ping, we used finite differences in space and two variable- 
order variable-step time-stepping methods: BDF (backward 
difference formulae) as a default solver and, the much slower, 
Adams method in those time intervals where there was no con- 
vergence achieved by the default solver. 

During the first few time-steps the front, originally at x = 
xq -c L, is quickly smoothed by surface tension, develops a 
capillary ridge and starts to move in the positive-;*: direction, 
thus the initial "dry patch" starts to grow. The concentra- 
tion field (j)(x,t) soon develops a peak at the dewetting front 
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that decays to the right. We need to resolve both the long 
decaying right hand tail of the capillary ridge h(x,t) during 
the convection-dominated stage and also the decay of the field 
(j)(x,t). Thus the length of the domain, L, should be larger than 
the decay lengths of h(x) and $(x,t), i.e. much larger than the 
lateral length scale a. On the other hand, we do not need to 
keep the dried region in our discretized domain where both 
h(x,t) and <j>(x,t) no longer change in time. Thus each time 
the front reaches some prescribed position x = x s hift> which 
is of the order of a, the simulation is automatically paused, 
and the computational domain is shifted (the x-coordinate also 
is shifted) to the right by the distance f s hif t = Jc s hift — x lst so 
that when the simulation is restarted the front becomes posi- 
tioned at a small distance x lst from the shifted boundary x = 0. 
We call this procedure shifting. The region through which the 
dewetting front travels up to the point jc s hift and the subsequent 
region where the capillary ridge needs to be resolved have the 
finest and equidistant grid-spacing. Then with increasing x 
the grid-spacing is gradually increased in order to obtain a 
more efficient numerical model. We always use the shift pe- 
riod /shift, that is an integer multiple of the finest grid-spacing 
so there is no interpolation involved in the region resolving the 
processes at the contact line and the capillary ridge. Several 
tens or hundreds of shifts usually happen during a typical sim- 
ulation that result in the deposition of 10-100 lines. Without 
shifting one would need to employ a much larger L and the 
region with finest grid-spacing would also need to be much 
larger than that we used. This would lead to a numerically 
intractable problem. 

After some performance and accuracy testing we deter- 
mined all the tuneable parameters of the numerics and then 
kept them fixed for all simulations published in this paper and 
also in the previous one.^ The values of the adjustable pa- 
rameters are: finest grid-spacing Ax = 8.5849, the period of 
shift f sh if t = 144Ajc, and the domain size L = 19230.2465. If 
we used the equidistant grid over the whole domain, the total 
number of grid-points would be 2241. However, because with 
our adapted grid the number of grid-points is 636. It could eas- 
ily be 10 times more than 2241 if we did not use the shifting 
procedure. 

A.2 Boundary conditions 

A semi-infinitely extended thin film evaporatively dewetting 
and depositing lines has a front that moves with an intrinsic 
velocity that can vary in time. The numerical solution method 
involves finite differences and so we are restricted to a finite 
domain. This does not cause a problem on the dried (left) 
side of the solution domain, where the unknown fields h(x,t) 
and hp(x,t) essentially do not evolve in time. Thus the left 
BCs at x = are simply symmetry conditions for h(x,t), i.e. 
d x h(x,t)\ x= Q = d xxx h(x,t)\ x= Q = 0; and for (j)(x,t) the fictitious 



additional grid point out of the domain is fixed at the value 
corresponding to the last value shifted out. These BCs are em- 
ployed because when the region at the left boundary is dried 
and (<p > 1), then h(x,t) = h\ and the concentration is jammed, 

<j)(x,t) = <j)(x). 

However, one needs to be more careful regarding the right 
BC at x = L. When we assume that the evaporation process is 
very slow we can treat the dynamics of the fields h(x,t) and 
(j)(x,t) as quasistationary, i.e. the front moves with a veloc- 
ity that varies very slowly in time. Then we can reduce our 
Eqs. ( 13 i and ( 14 1 to the approximate ODEs. After lineariz- 
ing about the constant film height h(x) = hi, we can solve the 
problem in the region far from the dewetting front to the right 
(wet) side and find stationary exponentially decaying profiles 
for h(x) and (j>(x). These profiles are used to prescribe the 
boundary conditions on the right side of our solution domain, 
x = L, which read 



d x h(x,t) = X(h(x,t) — hi), 
d xxx h(x,t)=X 3 (h(x,t)-h 2 ), 



(22) 
(23) 
(24) 



where <po is the initial concentration of the nanoparticles. Ini- 
tially, at the right boundary h(L,Q) = hi and <j>(L,Q) = <po and 
these boundary conditions reduce to setting all odd deriva- 
tives to zero. This is the same as in no-flux boundary con- 
ditions (see below Eqs. (|26]>). During the simulation always 
h(L,t) > hi and <j>(L,t) > 0o- The decay rate, X, is 



A = 



h\ 



<j>o[h(L,t) - h 2 ] 

/i 2 [<HM-<M 



(25) 



so that, apart from the constant parameters, the decay rate non- 
linearly depends on the values h(L,t) and Q(L,t). This means 
that the decay rate and the velocity of the moving front of a 
dilute suspension are different from the ones of a more dense 
suspension. 

Note, that the effect of an increasing average deposit thick- 
ness when decreasing the evaporation rate (as seen, in partic- 
ular, in the calculations at fixed 0o = 0.31 and varied Qq in 



Fig. 12 1 is due to a small systematic artefact connected to the 



simple approximations employed to derive the BCs (22 1— ( 24 1 
at the right boundary at x = L. These BCs correspond to a situ- 
ation with continuous systematic supply of additional solution 
during the convection-dominant regime and can be physically 
realised as some finite size effect. To further minimise this 
effect we used a large computational domain size L. 

Alternatively, we also performed several simulations with 
no-flux BCs: 

d x h(x, t) \ X=L = d xxx h(x,t)\ x=L = d x <j> (x,t)\ x=L = 0. (26) 

They prescribe zero flux through the right boundary and thus 
correspond to a different limit than the BCs that involve X. 
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Comparing results obtained with both types of right BCs 
shows no qualitative difference; see the comparison given in 
the next section. 

A.3 Remark on the robustness of our results 

In order to better utilise the computational resources, a num- 
ber of numerical tricks are employed: shifting, locally refined 
grid, and boundary conditions for decaying fields h(x,t) and 
(j>(x,t). Great care is taken to avoid any artefacts in the deposi- 
tion profiles due to these tricks. As a check, the computed pat- 
tern are in some cases verified by recomputing them with each 
of these refinements separately: (a) with different shift length, 
f s hift, that is further away from a possible resonance with 
the period of the pattern (180/144 of the original P s hift); (b) 
with refined grid, where in the region with equidistant finest 
grid-spacing we simply add one extra node between each two 
nodes of the original grid and accordingly we modify the rest 
of the grid without testing whether the grid-spacing outside 
the equidistant zone is optimal; (c) using the no-flux bound- 



ary conditions at the right boundary (see Appendix A. 2 1; (d) 
larger L (40/35 of the original L). For these tests, we se- 
lect three of the characteristic configurations used throughout 
this paper e.g. in Figs.[6] and [7] namely: the standard case 
(Qo,<h) = (4-64 x 10 _7 70.41), the case with significant neg- 
ative skewness (X2o,0o) = (4.64 x 10~ 6 ,0.41) and the large- 
amplitude long-period case (X2o,0o) = (4.64 x 10~ 7 ,0.31). It 
is found that the resulting patterns are qualitatively indepen- 
dent of these refinements, however we observe some differ- 
ences in the values of the computed measures, as follows: 

• The smallest errors are obtained for the test using a dif- 
ferent shift length: an error under 1.2% for the ampli- 
tude and period and under 2.6% for the relative width 
and skewness, except for the large-amplitude long-period 
case where the error in the period is 14%, in the relative 
width is 17.6% and in the skewness is 10.6%. 

• The errors from tests on the refined grid are: under 6% 
for the amplitude and under 15.8% for the period, except 
for the large-amplitude long-period case, where the error 
in the period is 55%. The errors for the relative width and 
the skewness are somewhat larger. 

• The remaining tests (using symmetric BC and larger L) 
give in some cases larger errors in the measures but 
changes in the BC or L may be thought of as varying 
the finite size effect. 

Thus, the numerical parameters: P s hifb the grid-spacing, BCs 
with A and L that were used for all presented simulations were 
fine tuned, so as to give the most reliable results possible, sub- 
ject to the computational resources available. 



B Measures for periodic patterns 

To study the changes in the morphology of the deposited pat- 
tern of periodic lines, /z p (x), we computed the following mea- 
sures in a selected spatial region containing the fully devel- 
oped periodic pattern: The basic measures are the amplitude, 
a, and the period, P, 



a = -E fl ' = -E( /i p(3'0-^p I ), 



n 



= 1 



^ n \ n 

^ =-Lfl = -L(^i-^), 



(27) 
(28) 



i=\ 



1=1 



where x; [y,] is the position of /th subsequent minimum [max- 
imum] of h p (x), respectively; hp. — j[/z p (x,) + /i p (x,- + i)] is the 
minimum height of the /th line (peak). By /th line we mean the 
portion of h p (x) that is defined on the interval The 
remaining measures are integral ones 



Xi 



Oi = 



h Pi (x)dx, 



hp j (x)xdx, 



73/ 



74/ 



hrt: = 




(x)(x — x,-) 2 dx, 



h Pi (x)(x—Xi) dx, 
h Pi (x)(x-Xi) 4 dx, 



(x)dx. 



(29) 
(30) 

(31) 
(32) 
(33) 

(34) 

(35) 



The given quantities for the /th peak are: A, is the area of the 
excess nanoparticle film, x, is the peak's centroid, a, is its stan- 
dard deviation, jq is its normalised width, 73; is its skewness, 
Y4i is its kurtosis, and /i P; . is its average height. The final mea- 
sures (that are used throughout this paper) hp, A, x, a, K, 73, 
74, and hp are obtained by averaging over n subsequent lines 
as is shown in the case of amplitude in Eq. (27 1. 

Note, that higher order moments such as a, 73, 74 give too 
much weight to the tail of a peak, especially in the case of 
large periods. For these, small numerical variations in the tail 
of a peak profile are greatly amplified in the value. Since our 
measures shall focus on the peak characteristics and not on the 
shallow valleys between them, we cut off the shallow part of 
the tails, i.e. we set 

hp (x) = max [/j p (x) — hp. — ca,-, 0] (36) 

where the cutoff parameter, c = 0.05, is used in computing the 
measures displayed in the main text. 
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